Spinpolarised CAP fails

[AleksBL]

cannot add imaginary part to dH-type hamiltonian when Hamiltonian object is spin-polarised. It works with is the Hamiltonian object the CAP is put into is not spin-polarised.

Additional snippet in the tbtncTools code from Gaetano:

    if spinpol == False:
        dH_CAP = si.Hamiltonian(geometry, dtype=np.complex128)
    else:
        dH_CAP = si.Hamiltonian(geometry, dtype=np.complex128, spin = si.Spin("polarized"))

call in script:

   dH = CAP(Hd3.geom, ['top', 'bottom'], dz_CAP = 20, spinpol = True)

error: TypeError Traceback (most recent call last) /tmp/ipykernel_23393/1798268003.py in <cell line: 71>() 69 E3.fois_gras(STM) 70 Device.fois_gras(Hd3) ---> 71 dH = CAP(Hd3.geom, ['top', 'bottom'], dz_CAP = 20, spinpol = True) 72 73 with sisl.get_sile(Device.dir+'/CAP.dH.nc', mode='w') as fh:

/work1/abalo/BC_TWIBIGR/Twist_Vacancy/tbtncTools.py in CAP(geometry, side, dz_CAP, write_xyz, zaxis, spinpol) 379 orbs = dH_CAP.geom.a2o(idx) # if you have just 1 orb per atom, then orb = ia 380 for orb,wz in zip(orbs, Wz): --> 381 dH_CAP[orb, orb] = complex(0, -wz) 382 CAP_list.append(idx) 383 #print(list2range_TBTblock(idx))

/dtu/sw/dcc/SL73/2021-nov/XeonE5-2660v3/gnu/11.2.0/python/packages/3.9.9/sisl/0.12.1/lib/python3.9/site-packages/sisl/sparse_geometry.py in setitem(self, key, val) 1298 key = tuple(key) + (dd,) 1299 self._def_dim = -1 -> 1300 self._csr[key] = val 1301 1302 @property

/dtu/sw/dcc/SL73/2021-nov/XeonE5-2660v3/gnu/11.2.0/python/packages/3.9.9/sisl/0.12.1/lib/python3.9/site-packages/sisl/sparse.py in setitem(self, key, data) 1163 # TODO we need some way to reduce these things 1164 # for integer stuff. -> 1165 data = asarray(data, self._D.dtype) 1166 # Places where there are nan will be set to zero 1167 data[isnan(data)] = 0

/dtu/sw/dcc/SL73/2021-nov/XeonE5-2660v3/gnu/11.2.0/python/packages/3.9.9/numba/0.54.1/lib/python3.9/site-packages/numpy/core/_asarray.py in asarray(a, dtype, order, like) 100 return _asarray_with_like(a, dtype=dtype, order=order, like=like) 101 --> 102 return array(a, dtype, copy=False, order=order) 103 104

TypeError: can't convert complex to float

3.9.9 (main, Mar 8 2022, 14:12:27) [GCC 11.2.0] 0.0.0.dev+4ce19bb38b609f4c584336a0cf6a180af05180a8

[zerothi]

This is because of a problem in the Spin class and that tano uses the dim variable.

Simply call Tano's routine as CAP(..., spin=2) which should work.

[AleksBL]

I didnt make your suggestion work, but im doing the following (obvious now) workaround:

Just writing the elements of one spin to another Hamiltonian object.

idx = [which_spin,2]
for I in range(3):
    for J in range(3):
        for i in range(no1):
            for j in range(no1):
                hd[i,j,(I,J,0)] = _hd[i,j,(I,J,0)][idx]
        for i in range(no2):
            for j in range(no2):
                hp[i,j,(I,J,0)] = _hp[i,j,(I,J,0)][idx]

This way is however horrendeously slow. Is there any faster way of doing this?

[zerothi]

Please show the entire code. I know tano made use of polarised cap.

[AleksBL]

import sisl
import numpy as np
from tbtncTools import CAP
import matplotlib.pyplot as plt
from siesta_python.funcs import SC
from siesta_python.siesta_python import SiP
from inspect import getsource
import os
print(os.environ['OMP_NUM_THREADS'])
print(getsource(CAP))
test  = True
dimx, dimy = 70,70
which_spin = 0

_hd   = sisl.get_sile('../Ham_from_Fei/def.TSHS').read_hamiltonian()
_hp   = sisl.get_sile('../Ham_from_Fei/pris.TSHS').read_hamiltonian()

_hd   = _hd.sub_orbital([i for i in range(_hd.na)], [3,6])
_hp   = _hp.sub_orbital([i for i in range(_hp.na)], [3,6])

hd    = sisl.Hamiltonian(_hd.geom, orthogonal = False)
hp    = sisl.Hamiltonian(_hp.geom, orthogonal = False)

#assert 1 == 0

no1,no2 = hd.no,hp.no
idx = [which_spin,2]
for _I in range(3):
    for _J in range(3):
        for i in range(no1):
            for j in range(no1):
                I = _I-1
                J = _J-1
                hd[i,j,(I,J,0)] = _hd[i,j,(I,J,0)][idx]
        for i in range(no2):
            for j in range(no2):
                hp[i,j,(I,J,0)] = _hp[i,j,(I,J,0)][idx]

for t1 in range(1000):
    if np.linalg.norm(t1*hp.cell[0])>dimx:
        break
for t2 in range(1000):
    if np.linalg.norm(t2*hp.cell[1])>dimy:
        break
print(t1,t2)

eps_tol = 0.32
idx_vac = [i for i in range(hp.na) if (np.linalg.norm(hp.xyz[i] - hd.xyz, axis = 1)>eps_tol).all()]
print(idx_vac)
assert len(idx_vac)==1
idx_vac = idx_vac[0]
r_vac   = hp.xyz[idx_vac]
Hp_tiled  = hp.tile(t2,1).tile(t1+2,0)
HEM = hp.copy()
HEP = hp.copy()
HEP.xyz += Hp_tiled.cell[0] - HEM.cell[0]

t1r,t2r    = 2,1
#IDX_REMOVE = hp.na * t2 * t1 + hp.na * (t2r -1) #+ idx_vac 
T          = hp.cell[0] * t1r + hp.cell[1] *t2r 
hd_r       = hd.copy()
hd_r.xyz  += T
dij        = np.linalg.norm(Hp_tiled.xyz - (r_vac + T), axis = 1)
IDX_REMOVE = np.where(dij == dij.min() )[0][0]
Hd         = Hp_tiled.remove([IDX_REMOVE])
idx        = [i for i in range(Hd.na) if (np.linalg.norm(hd_r.xyz - Hd.xyz[i],axis = 1)<eps_tol).any()]

Hd2 = Hd.replace(idx, hd_r, eps = eps_tol)
Hd2 = (Hd2 + Hd2.transpose())/2

C = Hd2.center() + np.array([5,0,0])
atom = sisl.Atom('H', orbitals=[sisl.Orbital(R = 1.1), sisl.Orbital(1.1)])

STM = sisl.Geometry([0, 0, 1], atoms=atom, sc=sisl.SuperCell([10, 10, 1], nsc=[1, 1, 3]))
dTj = np.linalg.norm(Hd2.xyz - C, axis = 1)
IDX = np.argsort(dTj)
for JJ in IDX:
    if Hd2.xyz[JJ,2]>3.0:
        idx = JJ
        break

STM = STM.move(Hd2.xyz[idx] - STM.xyz[0])
STM = sisl.Hamiltonian(STM,spin = sisl.Spin("polarized"))
#STM.construct([[0.1, 1.1], [0.0, -10.0]])
STM[0,0] = 0.0
STM[1,1] = 0.0
STM[0,0,(0,0,1)] = 10.0
STM[1,1,(0,0,1)] = 10.0
STM[0,0,(0,0,-1)] = 10.0
STM[1,1,(0,0,-1)] = 10.0

#Contour = np.linspace(1,-1, 3) + 1j * 1e-5

#assert 1 == 0

E1 = SiP(HEM.cell, HEM.xyz, HEM.toASE().numbers, pp_path = '../../pp',
         kp = [1,t2,1], semi_inf = '-a1',
         directory_name = 'E1', sl = 'E1',sm = 'E1',
         elec_bloch = [1,t2,1],
        )
E2 = SiP(HEP.cell, HEP.xyz, HEP.toASE().numbers, pp_path = '../../pp',
         kp = [1,t2,1], semi_inf = '+a1',
         directory_name = 'E2', sl = 'E2',sm = 'E2',
         elec_bloch = [1,t2,1],
        )
E3 = SiP(STM.cell, STM.xyz, STM.atoms.Z.copy(), pp_path = '../../pp',
        semi_inf = '+a3',
         directory_name = 'E3', sl = 'E3',sm = 'E3',
        )

Device = SiP(Hd2.cell.copy(), Hd2.xyz.copy(), Hd2.atoms.Z.copy(),
             Chem_Pot = [.0, .0, .0], 
             spin_pol = 'polarized',
             directory_name = 'Device',
             pp_path = '../../pp',
             kp = [1,1,1], 
             kp_tbtrans = [1,1,1],
             elecs = [E1,E2,E3],
             skip_distance_check = True,
             skip_spg_sp   =  True,
             trans_emin    = -1.,
             trans_emax    =  1.,
             trans_delta   =  1.,
             save_SE       =  True,
             save_SE_only  =  True,
             print_console =  True,
             #custom_tbtrans_contour = Contour,
            )
Device.find_elec_inds()
rear_idx = Device._rearange_indices
Hd3 = Hd2.sub(rear_idx)
E1.fois_gras(HEM)
E2.fois_gras(HEP)
E3.fois_gras(STM)
Device.fois_gras(Hd3)
dH = CAP(Hd3.geom, ['top', 'bottom'], dz_CAP = 20,#spin = 2
         #spinpol = True
        )

with sisl.get_sile(Device.dir+'/CAP.dH.nc', mode='w') as fh:
    fh.write_delta(dH)
#print(dH)

Device.fdf()
Device.write_more_fdf(['TBT.dH CAP.dH.nc'], name = 'TS_TBT')
Device.write_more_fdf(['TBT.Symmetry.TimeReversal false'], name = 'TS_TBT')
Device.mpi = ' '

Device.run_tbt_analyze_in_dir()
Device.run_tbtrans_in_dir(ORBITAL_CURRENT = True) # Fails here with some error code -40 in a netCDF file if the CAP is used

[zerothi]

I need the code that calls Tanos script, or even better could you make a small script that reproduces the problem? That would be easier to debug

[zerothi]

I don't understand that it does not work?

I did this:

import sisl as si
from tbtncTools import CAP

import sys

gr = si.geom.graphene() * (10, 10, 1)

dH = CAP(gr, ['top', 'bottom'], dz_CAP = 20, spin = 2)
print(dH)
print(dH._csr)
print(dH.dtype)

and it worked fine?

[AleksBL]

My version does not have this spin-keyword, I guess that is where the problem is.... xD

[zerothi]

I am using the latest from tanocalogero/mbare, I don't know if another one exists?

[zerothi]

This seemed to be related to incompatible versions. Will close issue.