Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fixed reading geometries from fdf + XV + STRUCT #779

Merged
merged 4 commits into from
May 29, 2024
Merged

fixed reading geometries from fdf + XV + STRUCT #779

merged 4 commits into from
May 29, 2024

Conversation

zerothi
Copy link
Owner

@zerothi zerothi commented May 29, 2024

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.

  • Closes Fix the basis when reading from XV #778
  • Added tests for new/changed functions?
  • Ran isort . and black . [24.2.0] at top-level
  • Documentation for functionality in docs/
  • Changes documented in CHANGELOG.md

Replaces #778

@@ -10,7 +10,7 @@
import pytest

import sisl
from sisl import Atom, Geometry, geom
from sisl import Atom, AtomUnknown, Geometry, geom

Check notice

Code scanning / CodeQL

Unused import Note test

Import of 'Atom' is not used.
Import of 'Geometry' is not used.
Import of 'AtomUnknown' is not used.
@tfrederiksen
Copy link
Contributor

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
Copy link
Owner Author

zerothi commented May 29, 2024

Could you share your files?

@zerothi
Copy link
Owner Author

zerothi commented May 29, 2024

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
Copy link
Owner Author

zerothi commented May 29, 2024

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

@tfrederiksen
Copy link
Contributor

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

zerothi added 2 commits May 29, 2024 13:41
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.

Signed-off-by: Nick Papior <[email protected]>
Copy link

codecov bot commented May 29, 2024

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 ⚠️
src/sisl/io/siesta/fdf.py 92.30% 1 Missing ⚠️
Additional details and impacted files
@@            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     

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

@tfrederiksen
Copy link
Contributor

69d0f3f works correctly for my system!

@pfebrer
Copy link
Contributor

pfebrer commented May 29, 2024

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
Copy link
Contributor

pfebrer commented May 29, 2024

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
Copy link
Owner Author

zerothi commented May 29, 2024

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

@zerothi
Copy link
Owner Author

zerothi commented May 29, 2024

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
Copy link
Contributor

pfebrer commented May 29, 2024

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 🤔

Signed-off-by: Nick Papior <[email protected]>
@zerothi
Copy link
Owner Author

zerothi commented May 29, 2024

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

@zerothi zerothi merged commit 2b299e3 into main May 29, 2024
3 checks passed
@zerothi
Copy link
Owner Author

zerothi commented May 29, 2024

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

@zerothi zerothi deleted the XV_fix_2nd branch May 29, 2024 18:53
@pfebrer
Copy link
Contributor

pfebrer commented May 29, 2024

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

@zerothi
Copy link
Owner Author

zerothi commented May 29, 2024

Omfg! Could you share it? :)

@zerothi
Copy link
Owner Author

zerothi commented May 29, 2024

Was it this? #779 (comment)

@zerothi
Copy link
Owner Author

zerothi commented May 29, 2024

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

@pfebrer
Copy link
Contributor

pfebrer commented May 29, 2024

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
Copy link
Owner Author

zerothi commented May 30, 2024

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

def test_fdf_multiple_simple(sisl_tmp):
.

@pfebrer
Copy link
Contributor

pfebrer commented May 30, 2024

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

@zerothi
Copy link
Owner Author

zerothi commented May 30, 2024

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
Copy link
Contributor

pfebrer commented May 30, 2024

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

@zerothi
Copy link
Owner Author

zerothi commented May 31, 2024

Should be fixed now! :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants