search

zerothi / sisl

Electronic structure Python package for post analysis and large scale tight-binding DFT/NEGF calculations
https://zerothi.github.io/sisl
Mozilla Public License 2.0
173 stars 57 forks source link

fixed reading geometries from fdf + XV + STRUCT #779

Closed zerothi closed 1 month ago

zerothi commented 1 month ago

The different geometries reading mechanisms where kind-of broken for various constalletations in the XV+fdf+struct files.

The species_as_Z was basicaly never going to works since there was no mechanism to counter it.

Instead the XV and STRUCT siles now handle basis from the perspective of being passed an atoms argument, from which it will fill in stuff based on data. If the atoms is a pure basis object (only correct species indices). Then it will replace without checking stuff. If it has the same length as the geometry it will check Z and whether indices are the same. This can potentially lead to errors if there is only 1 atom of each species, but on the other hand the species order should then handle the indices correctly.

There was also cases in fdf, XV and STRUCT reads that could potentially lead to reducing the basis. We should generally try to avoid this since there may be good reasons for having empty basis.

Replaces #778

tfrederiksen commented 1 month ago

This fails for me on the same files as reported in #778. With

g0 = sisl.get_sile("./RUN.fdf").read_geometry(output=True)
print(g0.atoms.atom)

I get this error message:


info:0: SislInfo: Replacing atom
  Atom{H, Z: 1, mass(au): 1.00794, maxR: -1.00000,
   Orbital{R: -1.00000, q0: 0.0}
  }
->
  Atom{H.opt88, Z: 1, mass(au): 1.01000, maxR: 2.85200,
   AtomicOrbital{1sZ1, q0: 1.0, SphericalOrbital{l: 0, R: 2.8179000000000167, q0: 1.0}},
   AtomicOrbital{1sZ2, q0: 0.0, SphericalOrbital{l: 0, R: 1.9918999999999998, q0: 0.0}},
   AtomicOrbital{2pyZ1P, q0: 0.0, SphericalOrbital{l: 1, R: 2.8520000000000043, q0: 0.0}},
   AtomicOrbital{2pzZ1P, q0: 0.0, SphericalOrbital{l: 1, R: 2.8520000000000043, q0: 0.0}},
   AtomicOrbital{2pxZ1P, q0: 0.0, SphericalOrbital{l: 1, R: 2.8520000000000043, q0: 0.0}}
  }
with a different number of orbitals!
Traceback (most recent call last):
  File "test.py", line 16, in <module>
    g0 = sisl.get_sile("./relax/RUN.fdf").read_geometry(output=True)
  File "/home/me/.local/lib/python3.10/site-packages/sisl/io/siesta/fdf.py", line 1443, in read_geometry
    v = getattr(self, f"_r_geometry_{f.lower()}")(*args, **kwargs)
  File "/home/me/.local/lib/python3.10/site-packages/sisl/io/siesta/fdf.py", line 1461, in _r_geometry_xv
    geom = xvSileSiesta(f).read_geometry(atoms=basis)
  File "/home/me/.local/lib/python3.10/site-packages/sisl/io/sile.py", line 707, in pre_open
    return func(self, *args, **kwargs)
  File "/home/me/.local/lib/python3.10/site-packages/sisl/messages.py", line 106, in wrapped
    return func(*args, **kwargs)
  File "/home/me/.local/lib/python3.10/site-packages/sisl/messages.py", line 106, in wrapped
    return func(*args, **kwargs)
  File "/home/me/.local/lib/python3.10/site-packages/sisl/messages.py", line 106, in wrapped
    return func(*args, **kwargs)
  File "/home/me/.local/lib/python3.10/site-packages/sisl/io/siesta/xv.py", line 160, in read_geometry
    _replace_basis(atms2, atoms)
  File "/home/me/.local/lib/python3.10/site-packages/sisl/io/siesta/_help.py", line 212, in _replace_basis
    if len(in_idx) > 0 and (np.allclose(in_idx, ref_idx) or only_basis):
  File "<__array_function__ internals>", line 200, in allclose
  File "/home/me/.local/lib/python3.10/site-packages/numpy/core/numeric.py", line 2270, in allclose
    res = all(isclose(a, b, rtol=rtol, atol=atol, equal_nan=equal_nan))
  File "<__array_function__ internals>", line 200, in isclose
  File "/home/me/.local/lib/python3.10/site-packages/numpy/core/numeric.py", line 2380, in isclose
    return within_tol(x, y, atol, rtol)
  File "/home/me/.local/lib/python3.10/site-packages/numpy/core/numeric.py", line 2361, in within_tol
    return less_equal(abs(x-y), atol + rtol * abs(y))
ValueError: operands could not be broadcast together with shapes (324,) (243,)
zerothi commented 1 month ago

Could you share your files?

zerothi commented 1 month ago

This fails for me on the same files as reported in #778. With

g0 = sisl.get_sile("./RUN.fdf").read_geometry(output=True)
print(g0.atoms.atom)

I get this error message:

info:0: SislInfo: Replacing atom
  Atom{H, Z: 1, mass(au): 1.00794, maxR: -1.00000,
   Orbital{R: -1.00000, q0: 0.0}
  }
->
  Atom{H.opt88, Z: 1, mass(au): 1.01000, maxR: 2.85200,
   AtomicOrbital{1sZ1, q0: 1.0, SphericalOrbital{l: 0, R: 2.8179000000000167, q0: 1.0}},
   AtomicOrbital{1sZ2, q0: 0.0, SphericalOrbital{l: 0, R: 1.9918999999999998, q0: 0.0}},
   AtomicOrbital{2pyZ1P, q0: 0.0, SphericalOrbital{l: 1, R: 2.8520000000000043, q0: 0.0}},
   AtomicOrbital{2pzZ1P, q0: 0.0, SphericalOrbital{l: 1, R: 2.8520000000000043, q0: 0.0}},
   AtomicOrbital{2pxZ1P, q0: 0.0, SphericalOrbital{l: 1, R: 2.8520000000000043, q0: 0.0}}
  }
with a different number of orbitals!
Traceback (most recent call last):
  File "test.py", line 16, in <module>
    g0 = sisl.get_sile("./relax/RUN.fdf").read_geometry(output=True)
  File "/home/me/.local/lib/python3.10/site-packages/sisl/io/siesta/fdf.py", line 1443, in read_geometry
    v = getattr(self, f"_r_geometry_{f.lower()}")(*args, **kwargs)
  File "/home/me/.local/lib/python3.10/site-packages/sisl/io/siesta/fdf.py", line 1461, in _r_geometry_xv
    geom = xvSileSiesta(f).read_geometry(atoms=basis)
  File "/home/me/.local/lib/python3.10/site-packages/sisl/io/sile.py", line 707, in pre_open
    return func(self, *args, **kwargs)
  File "/home/me/.local/lib/python3.10/site-packages/sisl/messages.py", line 106, in wrapped
    return func(*args, **kwargs)
  File "/home/me/.local/lib/python3.10/site-packages/sisl/messages.py", line 106, in wrapped
    return func(*args, **kwargs)
  File "/home/me/.local/lib/python3.10/site-packages/sisl/messages.py", line 106, in wrapped
    return func(*args, **kwargs)
  File "/home/me/.local/lib/python3.10/site-packages/sisl/io/siesta/xv.py", line 160, in read_geometry
    _replace_basis(atms2, atoms)
  File "/home/me/.local/lib/python3.10/site-packages/sisl/io/siesta/_help.py", line 212, in _replace_basis
    if len(in_idx) > 0 and (np.allclose(in_idx, ref_idx) or only_basis):
  File "<__array_function__ internals>", line 200, in allclose
  File "/home/me/.local/lib/python3.10/site-packages/numpy/core/numeric.py", line 2270, in allclose
    res = all(isclose(a, b, rtol=rtol, atol=atol, equal_nan=equal_nan))
  File "<__array_function__ internals>", line 200, in isclose
  File "/home/me/.local/lib/python3.10/site-packages/numpy/core/numeric.py", line 2380, in isclose
    return within_tol(x, y, atol, rtol)
  File "/home/me/.local/lib/python3.10/site-packages/numpy/core/numeric.py", line 2361, in within_tol
    return less_equal(abs(x-y), atol + rtol * abs(y))
ValueError: operands could not be broadcast together with shapes (324,) (243,)

Ok, this is actually one of the weird cases where you likely have two different H atoms, no? I should have a test with that one!

zerothi commented 1 month ago

I think I know how to amend, let me give it a try!

tfrederiksen commented 1 month ago

I have something like this:


%block ChemicalSpeciesLabel
  1    1   H.opt88
  2   79   Au.blyp
  3   79   Aus.blyp
  4    6   C.blyp
%endblock ChemicalSpeciesLabel

%block AtomicCoordinatesAndAtomicSpecies
x y z 1
x y z 2
x y z 3
x y z 4
x y z 1
%endblock AtomicCoordinatesAndAtomicSpecies
codecov[bot] commented 1 month ago

Codecov Report

Attention: Patch coverage is 98.80952% with 2 lines in your changes are missing coverage. Please review.

Project coverage is 87.36%. Comparing base (9ddceee) to head (760ee07).

Files Patch % Lines
src/sisl/io/_help.py 97.61% 1 Missing :warning:
src/sisl/io/siesta/fdf.py 92.30% 1 Missing :warning:
Additional details and impacted files ```diff @@ Coverage Diff @@ ## main #779 +/- ## ========================================== + Coverage 87.32% 87.36% +0.03% ========================================== Files 396 396 Lines 50592 50700 +108 ========================================== + Hits 44178 44292 +114 + Misses 6414 6408 -6 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

tfrederiksen commented 1 month ago

69d0f3f works correctly for my system!

pfebrer commented 1 month ago

Here this doesn't work:

XV_fail.zip

fdf:

NumberOfAtoms         2
NumberOfSpecies       2

%block ChemicalSpeciesLabel
  1   79  Au
  2    8   O
%endblock ChemicalSpeciesLabel

 LatticeConstant    1.000 Ang
 %block LatticeVectors
 10.475164246  0.000000000  0.000000000
 0.000000000   9.071764555  0.000000000
 0.000000000   0.000000000  44.467055785
 %endblock LatticeVectors
 AtomicCoordinatesFormat Ang
 %block AtomicCoordinatesAndAtomicSpecies
   0.0000000000   0.0000000000   0.0000000000   2
   2.6187910000   1.5119600000   0.0000000000   1
 %endblock AtomicCoordinatesAndAtomicSpecies

XV:

         19.795199425       0.000000000       0.000000000          0.000000000       0.000000000       0.000000000
          0.000000000      17.143157308       0.000000000          0.000000000       0.000000000       0.000000000
          0.000000000       0.000000000      84.030590492          0.000000000       0.000000000       0.000000000
         2
  2    8       0.000000000       0.000000000       0.000000000          0.000000000       0.000000000       0.000000000
  1    79       4.948799740       2.857191450       0.000000000          0.000000000       0.000000000       0.000000000

Which gives the following Z after reading with sisl:

python -c 'import sisl; print(sisl.get_sile("siesta.fdf").read_geometry(output=True).atoms.Z)'
[79  8]

For some reason it works if the basis is in a PAO.Basis block in the fdf file, but not if it is read from .ion.* files.

pfebrer commented 1 month ago

Even just reading the geometry from the fdf file (without any basis) is wrong:

import sisl

with open("geom.fdf", "w") as f:
    f.write("""
    NumberOfAtoms         2
    NumberOfSpecies       2

    %block ChemicalSpeciesLabel
      1   79  Au
      2    8   O
    %endblock ChemicalSpeciesLabel

     LatticeConstant    1.000 Ang
     %block LatticeVectors
     4  0  0
     0  10 0
     0  0  10
     %endblock LatticeVectors
     AtomicCoordinatesFormat Ang
     %block AtomicCoordinatesAndAtomicSpecies
       0 0 0  2
       2 0 0  1
     %endblock AtomicCoordinatesAndAtomicSpecies
     """)

sisl.get_sile("geom.fdf").read_geometry().atoms.Z
array([79,  8], dtype=int32)
zerothi commented 1 month ago

Thanks, my tests are obviously only testing species, not z... Damn... I'll have another stab at this!

zerothi commented 1 month ago

Actually, your test was exactly the kind of geometries that were the bottleneck, it would have worked if there were more atoms than species ;)

But, I will fix this!

pfebrer commented 1 month ago

Hmm I have another one that fails and has 260 atoms and 5 species. I thought this was a minimal example of it, but let's see 😅

Actually if you add another oxygen to the example that I provided it is also wrong 🤔

zerothi commented 1 month ago

I added your Au, C example as a test, it passed, could you try the latest commit out?

zerothi commented 1 month ago

I'll just merge, it seems quite important ;) @pfebrer let me know if there are still issues!

pfebrer commented 1 month ago

This example (with the ion files) still doesn't work :sweat_smile:

zerothi commented 1 month ago

Omfg! Could you share it? :)

zerothi commented 1 month ago

Was it this? https://github.com/zerothi/sisl/pull/779#issuecomment-2137678736

zerothi commented 1 month ago

I'll see if I get time tomorrow... Sorry!

pfebrer commented 1 month ago

Yes sorry, forgot to link it. There the species are still reversed.

And if I add one extra oxygen I get:

warn:0: SislWarning: Trying to replace atom <sisl.Atom O, Z=8, M=16.0, maxR=3.1449000000000105, no=13> with <sisl.Atom O, Z=8, M=16.0, maxR=3.1449000000000105, no=13>, but they don't share the same atoms in the geometry.
warn:0: SislWarning: Trying to replace atom <sisl.Atom Au, Z=79, M=196.97, maxR=3.3210000000000024, no=15> with <sisl.Atom Au, Z=79, M=196.97, maxR=3.3210000000000024, no=15>, but they don't share the same atoms in the geometry.

I don't understand why this simple case (reading XV from a fdfSileSiesta) should substitute the basis by checking things. The fdf contains a dictionary (or list) of species:

%block ChemicalSpeciesLabel
     1 79 Au
     2 8 O
%endblock

and the first column of the XV refers explicitly to the key/index of the species as defined in ChemicalSpeciesLabels. So there's nothing to check there, no? We just need to use whatever species the XV file says...

zerothi commented 1 month ago

Yes sorry, forgot to link it. There the species are still reversed.

And if I add one extra oxygen I get:

warn:0: SislWarning: Trying to replace atom <sisl.Atom O, Z=8, M=16.0, maxR=3.1449000000000105, no=13> with <sisl.Atom O, Z=8, M=16.0, maxR=3.1449000000000105, no=13>, but they don't share the same atoms in the geometry.
warn:0: SislWarning: Trying to replace atom <sisl.Atom Au, Z=79, M=196.97, maxR=3.3210000000000024, no=15> with <sisl.Atom Au, Z=79, M=196.97, maxR=3.3210000000000024, no=15>, but they don't share the same atoms in the geometry.

I don't understand why this simple case (reading XV from a fdfSileSiesta) should substitute the basis by checking things. The fdf contains a dictionary (or list) of species:

%block ChemicalSpeciesLabel
     1 79 Au
     2 8 O
%endblock

and the first column of the XV refers explicitly to the key/index of the species as defined in ChemicalSpeciesLabels. So there's nothing to check there, no? We just need to use whatever species the XV file says...

But it seems you have some extra files, I would guess it actually reads the basis from somewhere else, I added the test for your case, and it passes fine, see https://github.com/zerothi/sisl/blob/58f4758438a66f598e6140826a87571c3e047c12/src/sisl/io/siesta/tests/test_fdf.py#L686.

pfebrer commented 1 month ago

Yes, https://github.com/zerothi/sisl/pull/779#issuecomment-2137678736 contains the XV and the ion files (the test is only the FDF)

zerothi commented 1 month ago

Yes, #779 (comment) contains the XV and the ion files (the test is only the FDF)

Thanks! I had missed that ;)

It should be fixed in the latest commit! A problem only occurring when reading basis from ion files... :sigh:

pfebrer commented 1 month ago

Btw, now the tests are full of warnings about replacing atoms 😅

zerothi commented 1 month ago

Should be fixed now! :)