translate2uc for sparse matrices with associated geometries.

[pfebrer]

As discussed in #583, this PR creates a function that translates all atoms of a sparse matrix to the unit cell, updating both the coordinates and the indices of the matrix elements.

The current state of the code contains only the simplest possible behavior, implemented in the most straightforward way I could think of. And only for SparseOrbital, not yet for SparseAtom. This is so that you can easily give it a look and see if there is some conceptual mistake. Also, I am aware that this is not the most time and memory efficient implementation possible, but I can optimize it later (suggestions accepted :)). Also I don't know what is the approppiate ratio optimization/readability for this.

I thought about it, and I see no way of using the translate_columns without making the code more complex and unreadable, because the translation of each matrix element depends on both the columns and rows. So, a column translates differently depending on which row are we.

Demo

I did a very simple test with a hydrogen molecule, which I generated like:

import sisl
from ase.build import molecule
import numpy as np
import plotly.express as px

geom = sisl.Geometry.new(molecule("H2"))

geom.cell[0,0] = 10
geom.cell[1,1] = 10
geom.cell[2,2] = 10

geom.xyz += 5

Then I ran a SIESTA calculation for this molecule as is and the same molecule with the first atom displaced one cell in X. I used ForceAuxCell t and defaults for the rest of inputs.

I read the DMs and plot them:

uc_dm = sisl.get_sile("original/RUN.fdf").read_density_matrix()
sc_dm = sisl.get_sile("displaced/RUN.fdf").read_density_matrix()

px.imshow(uc_dm.tocsr().toarray()).show()
px.imshow(sc_dm.tocsr().toarray()).show()

Then, I translate the second DM to the unit cell:

moved_dm = sc_dm.translate2uc()

px.imshow(moved_dm.tocsr().toarray()).show()

All the heavy lifting is done by the translate_atoms_sc method, which receives the cell translations of all atoms. So you can displace the atoms differently and it also works. Here I displace the atom that was in the unit cell to the same cell as the displaced atom:

moved_dm = sc_dm.translate_atoms_sc([[0, 0, 0], [1,0,0]])

moved_dm.plot(axes="xz").show()
px.imshow(moved_dm.tocsr().toarray()).show()

[pfebrer]

One thing that I'm not sure how to do is to re-sort the column indices (that is, I don't know if there is already a method that does that). Because I modified them but I left them unsorted. scipy doesn't have a problem with it, but probably is better to sort them?

[zerothi]

Lets start with the first method, and work from there.

[zerothi]

Probably the interface should be equivalent to the Geometry.translate2uc, it would also allow selection of atoms.

[pfebrer]

Yes 👍 This was just to make it as simple as possible for now so that the code was easier to review

[zerothi]

Hmm... I don't particularly like this method. Couldn't we start with the translate2uc with the atoms specification, and if there is a need for the translate_atoms_sc then reconsider?

I see this method as being the equivalent of a translate method, since the geometry_coords can change the coords at will.

[pfebrer]

I will remove the geometry_coords argument, it was just to avoid recomputing the coordinates in the particular case of translate2uc but now I see that it's exactly the same amount of computation whether I precompute or compute inside the method.

I can also make the method private for now. But I thought it was a nice separation of functionality since the only specificity of translate2uc is to know how much unit cells is each atom translated. And since the translation of the matrix elements is not trivial to implement I figured it would be nice to keep it like this.

[pfebrer]

See now, you still don't like translate_atoms_sc? Not the name, but the concept 😅

[codecov]

Codecov Report

Merging #585 (9862512) into main (adb10c5) will increase coverage by 0.05%.
The diff coverage is 95.22%.

@@            Coverage Diff             @@
##             main     #585      +/-   ##
==========================================
+ Coverage   87.28%   87.33%   +0.05%     
==========================================
  Files         375      373       -2     
  Lines       47564    47676     +112     
==========================================
+ Hits        41518    41640     +122     
+ Misses       6046     6036      -10     

Impacted Files	Coverage Δ
src/sisl/geometry.py	86.14% <25.00%> (+0.25%)	⬆️
src/sisl/io/_multiple.py	88.37% <88.37%> (+1.55%)	⬆️
src/sisl/sparse_geometry.py	90.37% <90.69%> (+0.20%)	⬆️
src/sisl/io/sile.py	84.82% <96.15%> (ø)
src/sisl/io/orca/stdout.py	97.60% <99.12%> (-0.01%)	⬇️
src/sisl/io/openmx/md.py	100.00% <100.00%> (ø)
src/sisl/io/orca/tests/test_stdout.py	100.00% <100.00%> (ø)
src/sisl/io/siesta/ani.py	100.00% <100.00%> (ø)
src/sisl/io/siesta/tests/test_ani.py	100.00% <100.00%> (ø)
src/sisl/io/tests/test_xsf.py	100.00% <100.00%> (ø)
... and 5 more

... and 7 files with indirect coverage changes

[pfebrer]

I just checked that doing a dm.translate2uc before computing the density with the current algorithm fixes cases that returned a wrong density :) :)

Also I checked that the translation (even with the current perhaps naive implementation) is for free compared to the density computation. In my test system 35.5 ms versus 1.5min.

[pfebrer]

This is now fully functional, with tests that cover everything :)

After you review I can squash everything into a single commit.

[zerothi]

@pfebrer could you have a look? I would think this is fine (made the translate_atoms_sc a private function).

[zerothi]

Also, once in you should check that this works faster due to cb7e3cb

[zerothi]

Hold on, something is missing, a test for sparse orbital

[pfebrer]

Looks good to me 👍