Computing DM from H through eigenstates.

[pfebrer]

Trying to close #58.

I don't know exactly where do you want to put this functionality, but I open this PR with a working implementation that can be then reorganized as you wish whenever you want.

Basically, I introduced a compute_dm function that from a BrillouinZone object with an associated hamiltonian computes the DM. You can provide a custom distribution function to weight the eigenstates as I show below.

First, here's a test with a bunch of calculations of a Pt-Au chain with all spin variants at gamma and with 3 k points (files):

from pathlib import Path
import numpy as np

import sisl
from sisl.physics.compute_dm import compute_dm

# Define the root of the tests
root = Path("DM_tests")

results = []
# Loop over all directories in the tests directory
for test_dir in root.glob("*"):
    # Get fdf
    fdf = sisl.get_sile(test_dir / "RUN.fdf")
    # Read hamiltonian
    H = fdf.read_hamiltonian()
    
    # Read brillouin zone so that we can use exactly the same for our calculation
    kp = sisl.get_sile(fdf.file.parent / (fdf.get("SystemLabel", "siesta") + ".KP"))
    bz = kp.read_brillouinzone(H)
    
    # Read siesta's DM 
    siesta_DM = fdf.read_density_matrix()
    # Compute DM in sisl
    DM = compute_dm(bz, eta=False)
    
    # Calculate the difference
    errors = abs(siesta_DM - DM)
    
    # Add the maximum difference to the list of results
    results.append(
        (test_dir.name, np.max(errors._csr.data[:, :]))
    )
    
results

Which gives the following maximum differences:

[('unpol_k', 3.2786941075446663e-06),
 ('pol_gamma', 1.3595180536896123e-11),
 ('nc_gamma', 1.3598899784028617e-11),
 ('unpol_gamma', 2.398574228124062e-09),
 ('pol_k', 1.0009657427922924e-06),
 ('nc_k', 1.2885270135321036e-06),
 ('so_k', 1.1579284937557333e-06),
 ('so_gamma', 1.523530190894462e-10)]

An example of using it to get the LDOS:

root = Path("DM_tests")
fdf = sisl.get_sile(root / "pol_gamma" / "RUN.fdf")

# Read Hamiltonian and Brillouin zone
H = fdf.read_hamiltonian()
kp = sisl.get_sile(fdf.file.parent / (fdf.get("SystemLabel", "siesta") + ".KP"))
bz = kp.read_brillouinzone(H)

# Define a custom function to select the states that we wish
def in_range(x):
    return np.int32((x > -0.6) & (x < -0.4)).astype(float)

# Compute the DM for those states.
DM_LDOS = compute_dm(bz, occ_distribution=in_range, eta=False)

# Project the density on a grid
ldos_grid = sisl.Grid(0.1,geometry=DM.geometry)
DM_LDOS.density(ldos_grid)

# And plot it
ldos_grid.plot(axes="yx", zsmooth="best", nsc=[2, 2, 1])

Which is basically the same as the LDOS generated by SIESTA selecting the same states:

siesta_ldos = sisl.get_sile(fdf.file.parent / (fdf.get("SystemLabel", "siesta") + ".LDOS")).read_grid()
siesta_ldos.plot(axes="yx", zsmooth="best", nsc=[2,2,1])

Cheers!

[pfebrer]

By the way, I think it requires Cython 3.0 to compile properly as it is now, because of the cython.fused_type call.

[pfebrer]

This works (It's very easy for you to reproduce the test above and you don't need to compile with cython to check that it works since it's in python syntax). So whenever you want you can take this and reorganize it as you wish :) It would be very nice to have in sisl to finally allow LDOS(E) maps in combination with #496!