Inconsistency: the eigenvalues generated by Openmx code and TB2J-OpenMX
Closed this issue · 2 comments
When I used "TB2J-OpenMX/TB2J_OpenMX/ffiparser.py" to obtain the eigenvalues of a given k-point for the case of spin-orbit coupling calculations, I found that the eigenvalues obtained by TB2J-Openmx do not seem to be exactly equal to those obtained by Openmx code. TB2J-Openmx can work well for both the collinear and noncollinear cases.
For the example of "SrMnO3_FM_SOC":
k1= -0.41667 k2= -0.41667 k3= -0.41666
#band_index openmx_output TB2J-Opemx
1 -3.36190765415304 -3.361907667786186
2 -3.23140562557682 -3.231405644925162
3 -2.22902420861756 -2.228890673650604
4 -2.21444725079230 -2.215191732896975
5 -2.19501901702534 -2.195019050505414
6 -2.09419717272430 -2.094376973121326
7 -2.07974707772342 -2.080846583352733
8 -2.06773924016769 -2.067739250912386
9 -1.53490992796813 -1.534910354359798
10 -1.53464457466509 -1.534644719311316
11 -0.97657324061658 -0.977773083749660
12 -0.97415093549122 -0.974392952696173
13 -0.96902940696838 -0.969455363568774
14 -0.96853085804130 -0.968614391418724
15 -0.96642637880680 -0.966618152983887
16 -0.96588863049486 -0.966060411550144
17 -0.91192682008785 -0.925599823965553
18 -0.91154921990118 -0.924917820854717
19 -0.87583855698854 -0.890634411180473
20 -0.87512028165428 -0.885611349751454
21 -0.87216255246841 -0.874020166374503
22 -0.87148260939178 -0.872805132913553
23 -0.58248102861994 -0.582529609209740
24 -0.57563853391455 -0.575644305987460
25 -0.56740878564631 -0.570055447126084
26 -0.56560536390270 -0.565922132567649
27 -0.56500822133814 -0.565018833145592
28 -0.56378276404950 -0.564487765099754
29 -0.56325970868256 -0.563832455528609
30 -0.52270441687692 -0.522703757053155
31 -0.52265949977911 -0.522654952626849
32 -0.48889847351589 -0.491429663290982
33 -0.48316205574347 -0.485639174084016
34 -0.48159644500647 -0.482000673327639
35 -0.38089302414923 -0.382656872328660
36 -0.37771892333616 -0.379260351037036
37 -0.37744245771822 -0.377653986158365
38 -0.36333638693531 -0.363369712197737
39 -0.36302197627914 -0.363034812622270
40 -0.35895588583096 -0.358936267702248
41 -0.35740621604264 -0.357444556794187
42 -0.35714901502986 -0.357158401309471
43 -0.35296626535397 -0.352964052843812
44 -0.22484068040839 -0.225629642034248
45 -0.22374301434962 -0.224685879112824
46 -0.22289286512817 -0.222907250893914
47 -0.20096402211111 -0.200977920468490
The reason for the above inconsistency has been discovered. This is because the way employed in TB2J_Openmx (ffiparser.py) to construct the H and S matrix is different from the one in Openmx.
In ffiparser.py:
def parse_scfoutput(self):
argv0 = ffi.new("char[]", b"")
argv = ffi.new("char[]", bytes(self.fname, encoding='ascii'))
lib.read_scfout([argv0, argv])
lib.prepare_HSR()
self.ncell = lib.TCpyCell + 1
self.natom = lib.atomnum
self.norbs = np.copy(
asarray(ffi, lib.Total_NumOrbs, self.natom + 1)[1:])
if lib.SpinP_switch == 3:
self.non_collinear = True
elif lib.SpinP_switch == 1:
self.non_collinear = False
else:
raise ValueError(
" The value of SpinP_switch is %s. Can only get J from collinear and non-collinear mode."%lib.SpinP_switch)
fnan = asarray(ffi, lib.FNAN, self.natom + 1)
natn = []
for iatom in range(self.natom):
natn.append(asarray(ffi, lib.natn[iatom + 1], fnan[iatom + 1] + 1))
ncn = []
for iatom in range(self.natom):
ncn.append(asarray(ffi, lib.ncn[iatom + 1], fnan[iatom + 1] + 1))
# atv
# x,y,and z-components of translation vector of
#periodically copied cells
#size: atv[TCpyCell+1][4];
atv = []
for icell in range(self.ncell):
atv.append(asarray(ffi, lib.atv[icell], 4))
atv = np.copy(np.array(atv))
atv_ijk = []
for icell in range(self.ncell):
atv_ijk.append(asarray(ffi, lib.atv_ijk[icell], 4))
atv_ijk = np.array(atv_ijk)
self.R = atv_ijk[:, 1:]
tv = []
for i in range(4):
tv.append(asarray(ffi, lib.tv[i], 4))
tv = np.array(tv)
self.cell = np.copy(tv[1:, 1:]) * Bohr
rtv = []
for i in range(4):
rtv.append(asarray(ffi, lib.rtv[i], 4))
rtv = np.array(rtv)
self.rcell = np.copy(rtv[1:, 1:])
Gxyz = []
for iatom in range(self.natom):
Gxyz.append(asarray(ffi, lib.Gxyz[iatom + 1], 60))
self.positions = np.copy(np.array(Gxyz)[:, 1:4]) * Bohr
self.MP = np.copy(asarray(ffi, lib.MP, self.natom + 1)[1:])
self.norb = lib.T_NumOrbs
norb = self.norb
if self.non_collinear:
HR = np.zeros([self.ncell, 4, lib.T_NumOrbs, lib.T_NumOrbs])
for iR in range(0, self.ncell):
for ispin in range(lib.SpinP_switch + 1):
for iorb in range(lib.T_NumOrbs):
HR[iR, ispin,
iorb, :] = asarray(ffi, lib.HR[iR][ispin][iorb],
norb)
HR_imag = np.zeros([self.ncell, 3, lib.T_NumOrbs, lib.T_NumOrbs])
for iR in range(0, self.ncell):
for ispin in range(3):
for iorb in range(lib.T_NumOrbs):
HR_imag[iR, ispin, iorb, :] = asarray(
ffi, lib.HR_imag[iR][ispin][iorb], norb)
self.H = np.zeros(
[self.ncell, lib.T_NumOrbs * 2, lib.T_NumOrbs * 2],
dtype=complex)
for iR in range(self.ncell):
"""
# up up
self.H[iR, ::2, ::2] = HR[iR, 0, :, :] + 1j * HR_imag[iR, 0, :, :]
# up down
self.H[iR, ::2,
1::2] = HR[iR, 2, :, :] + 1j * (HR[iR, 3, :, :] +
HR_imag[iR, 2, :, :])
# down up
self.H[iR,
1::2, ::2] = HR[iR, 2, :, :] - 1j * (HR[iR, 3, :, :] +
HR_imag[iR, 2, :, :])
# down down
self.H[iR, 1::2, 1::2] = HR[iR, 1, :, :] + 1j * HR_imag[iR, 1, :, :]
"""
#following the way in openmx code:
#up up
self.H[iR, 0:norb, 0:norb] = HR[iR, 0, :, :] + 1j * HR_imag[iR, 0, :, :]
#up down
self.H[iR, 0:norb,
norb:2*norb] = HR[iR, 2, :, :] + 1j * (HR[iR, 3, :, :] +
HR_imag[iR, 2, :, :])
# down up
self.H[iR,
norb:2*norb, 0:norb] = HR[iR, 2, :, :] - 1j * (HR[iR, 3, :, :] +
HR_imag[iR, 2, :, :])
# down down
self.H[iR, norb:2*norb, norb:2*norb] = HR[iR, 1, :, :] + 1j * HR_imag[iR, 1, :, :]
"""
SR = np.zeros([self.ncell, lib.T_NumOrbs, lib.T_NumOrbs])
for iR in range(0, self.ncell):
for iorb in range(lib.T_NumOrbs):
SR[iR, iorb, :] = asarray(ffi, lib.SR[iR][iorb], norb)
self.S = np.kron(np.eye(2), SR)
"""
else: # collinear
HR = np.zeros([self.ncell, 4, lib.T_NumOrbs, lib.T_NumOrbs])
for iR in range(0, self.ncell):
for ispin in range(lib.SpinP_switch + 1):
for iorb in range(lib.T_NumOrbs):
HR[iR, ispin,
iorb, :] = asarray(ffi, lib.HR[iR][ispin][iorb],
norb)
self.H = np.zeros(
[self.ncell, lib.T_NumOrbs * 2, lib.T_NumOrbs * 2],
dtype=complex)
"""
for iR in range(self.ncell):
# up up
self.H[iR, ::2, ::2] = HR[iR, 0, :, :]
# down down
self.H[iR, 1::2, 1::2] = HR[iR, 1, :, :]
"""
for iR in range(self.ncell):
# up up
self.H[iR, 0:norb, 0:norb] = HR[iR, 0, :, :]
# up down
self.H[iR, 0:norb, norb:2*norb] = HR[iR, 2, :, :]
# down up
self.H[iR, norb:2*norb, 0:norb] = HR[iR, 2, :, :]
# down down
self.H[iR, norb:2*norb, norb:2*norb] = HR[iR, 1, :, :]
self.efermi = lib.ChemP * Ha
self.H *= Ha
SR = np.zeros([self.ncell, lib.T_NumOrbs, lib.T_NumOrbs])
for iR in range(0, self.ncell):
for iorb in range(lib.T_NumOrbs):
SR[iR, iorb, :] = asarray(ffi, lib.SR[iR][iorb], norb)
#self.S = np.kron(SR, np.eye(2))
self.S = np.kron(np.eye(2), SR)
lib.free_HSR()
lib.free_scfout()
print("Loading from scfout file OK!")
This has been discussed in
https://groups.google.com/g/tb2j/c/hAhJx0-kc48.
The two implementations are equivalent but leads to different numerical noise. The noise is from the non-Hermitian Hamiltonian in the OpenMX Calculation, which might need to be improved by using e.g. more k-points or more grid points.