Read/write of TSDE files

[tfrederiksen]

Describe the issue I wish to transform a .TSDE file, obtained from an unpolarized transiesta calculation, into a new one to be used in a spin-polarized run. However, I face some problems with the following script for a simple 1D test system:

1. Transformation spin-polarized to spin-polarized is OK As a sanity check, the following read-write script works well if one reads spin-polarized transiesta outputs (but not what I want):

import sisl

fh = sisl.get_sile('./RUN.fdf')
ef = fh.read_fermi_level()
DM = fh.read_density_matrix()
EDM = fh.read_energy_density_matrix()

DM = DM.transform(spin='polarized') # unnecessary in this case, because input is from a polarized run
EDM = EDM.transform(spin='polarized')

fh2 = sisl.get_sile('device.DM')
fh2.write_density_matrix(DM)
fh3 = sisl.get_sile('device.TSDE')
fh3.write_density_matrices(DM, EDM, Ef=ef)

I.e., read and overwrite of the .DM and .TSDE files makes a good starting point for another transiesta run:

$ more RUN.out
...
        iscf     Eharris(eV)        E_KS(eV)     FreeEng(eV)     dDmax    Ef(eV) dHmax(eV)
ts-scf:    1   -22318.455461   -22318.455562   -22318.455562  0.000014 -5.319343  0.001322
     spin moment: S , {S} =    0.00000       0.0       0.0   0.00000
timer: Routine,Calls,Time,% = TS             1       0.244   0.84
ts-q:         D        E1        C1        E2        C2        dQ   Qup-Qdn
ts-q:    53.592    43.567     0.410    43.566     0.417 0.438E-02 0.000E+00
ts-Vha:  0.10586E+00 eV
ts-scf:    2   -22318.457918   -22318.456729   -22318.456729  0.000048 -5.319343  0.004031
     spin moment: S , {S} =    0.00000       0.0       0.0   0.00000
ts-q:         D        E1        C1        E2        C2        dQ   Qup-Qdn
ts-q:    53.592    43.567     0.410    43.566     0.417 0.412E-02 0.000E+00
ts-Vha:  0.10709E+00 eV
ts-scf:    3   -22318.453974   -22318.455366   -22318.455366  0.000035 -5.319343  0.001118
     spin moment: S , {S} =    0.00000       0.0       0.0   0.00000
ts-q:         D        E1        C1        E2        C2        dQ   Qup-Qdn
ts-q:    53.592    43.567     0.410    43.566     0.417 0.419E-02 0.000E+00
ts-Vha:  0.10679E+00 eV
ts-scf:    4   -22318.456054   -22318.455706   -22318.455706  0.000008 -5.319343  0.000070
     spin moment: S , {S} =    0.00000       0.0       0.0   0.00000

2. Transformation unpolarized to spin-polarized fails Now, if I instead change the filehandle to look up the outputs from a converged, unpolarized run, i.e., changing the second line to

fh = sisl.get_sile('../TSrun_unpolarized/RUN.fdf')

which is basically what I wanted, the charges appear to be completely wrong with the .DM and .TSDE written by sisl:

$ more RUN.out
...
        iscf     Eharris(eV)        E_KS(eV)     FreeEng(eV)     dDmax    Ef(eV) dHmax(eV)
ts-scf:    1    60436.908351    -3832.241013    -3832.241013  2.000150 -5.319351939.969334
     spin moment: S , {S} =    0.00000       0.0       0.0   0.00000
timer: Routine,Calls,Time,% = TS             1       0.217   0.81
ts-q:         D        E1        C1        E2        C2        dQ   Qup-Qdn
ts-q:    17.353    87.134    -0.843    87.131    -0.523 0.492E+02 0.000E+00
ts-Vha:  0.42895E+03 eV
ts-scf:    2    13956.640735     1388.102402     1388.102402  2.508486 -5.319351**********
     spin moment: S , {S} =    0.00000       0.0       0.0   0.00000
ts-q:         D        E1        C1        E2        C2        dQ   Qup-Qdn
ts-q:     0.759    87.134    -1.062    87.131    -0.879 0.320E+02 0.000E+00
ts-Vha:  0.50424E+03 eV
ts-scf:    3    25851.044516     8597.095819     8597.095819  2.500046 -5.319351**********
     spin moment: S , {S} =    0.00000       0.0       0.0   0.00000
ts-q:         D        E1        C1        E2        C2        dQ   Qup-Qdn
ts-q:     0.005    87.134    -0.006    87.131    -0.007 0.332E+02 0.000E+00
ts-Vha:  0.50204E+03 eV
ts-scf:    4     7916.994467     8550.825876     8550.825876  0.362378 -5.319351770.460815
     spin moment: S , {S} =    0.00000       0.0       0.0   0.00000
ts-q:         D        E1        C1        E2        C2        dQ   Qup-Qdn
ts-q:     0.003    87.134    -0.004    87.131    -0.005 0.332E+02 0.000E+00
ts-Vha:  0.50204E+03 eV
ts-scf:    5     8549.494567     8550.938171     8550.938171  0.004230 -5.319351400.011045
     spin moment: S , {S} =    0.00000       0.0       0.0   0.00000
ts-q:         D        E1        C1        E2        C2        dQ   Qup-Qdn
ts-q:     0.003    87.134    -0.004    87.131    -0.005 0.332E+02 0.000E+00
ts-Vha:  0.50204E+03 eV
ts-scf:    6     8550.884241     8550.935553     8550.935553  0.000282 -5.319351349.250192
     spin moment: S , {S} =    0.00000       0.0       0.0   0.00000

3. Reference unpolarized output For reference, this is the convergence loop from the unpolarized run:

        iscf     Eharris(eV)        E_KS(eV)     FreeEng(eV)     dDmax    Ef(eV) dHmax(eV)
ts-scf:    1   -22318.453633   -22318.454524   -22318.454524  0.000054 -5.319351  0.002606
timer: Routine,Calls,Time,% = TS             1       0.107   0.44
ts-q:         D        E1        C1        E2        C2        dQ
ts-q:    53.592    43.567     0.410    43.566     0.417 0.486E-02
ts-Vha:  0.10363E+00 eV
ts-scf:    2   -22318.463456   -22318.458947   -22318.458947  0.000185 -5.319351  0.010229
ts-q:         D        E1        C1        E2        C2        dQ
ts-q:    53.592    43.567     0.410    43.566     0.417 0.418E-02
ts-Vha:  0.10681E+00 eV
ts-scf:    3   -22318.451810   -22318.455416   -22318.455416  0.000149 -5.319351  0.000092
ts-q:         D        E1        C1        E2        C2        dQ
ts-q:    53.592    43.567     0.410    43.566     0.417 0.418E-02
ts-Vha:  0.10682E+00 eV
ts-scf:    4   -22318.455383   -22318.455399   -22318.455399  0.000002 -5.319351  0.000088

Any clues what could be the issue here?

Version details

>>> import sys
>>> print(sys.version)
3.10.6 (main, Mar 10 2023, 10:55:28) [GCC 11.3.0]
>>> import sisl
>>> print(sisl.__version__)
0.13.1.dev119+g85435d44

(Siesta Version : 4.1.5)

[zerothi]

Yes, when upscaling a density matrix, you need to halve the density, since the unpolardzed to polarized just copies the elements. So

DM = DM.transform(...) / 2

Should work, also note the charges in the electrode which should hint at this.

[tfrederiksen]

You're absolutely correct. With the scaling the script works, thanks!

[tfrederiksen]

For completeness, the EDM should in principle also be divided by a factor two, right? (In my test case such a factor only seems to affect the calculated forces in the electrodes).

[zerothi]

Yes, for density related stuff it is important, for the Hamiltonian it is not. :)

[zerothi]

For completeness, the EDM should in principle also be divided by a factor two, right? (In my test case such a factor only seems to affect the calculated forces in the electrodes).

Probably because of having DM.update cross which means that the EDM in the device gets updated. :)