diff --git a/CHANGELOG.md b/CHANGELOG.md index 1a4826aa3c..66f2bfab03 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,13 +24,12 @@ we hit release version 1.0.0. - Siesta sparse matrices could in some cases set wrong diagonal components - orbital quantum numbers from HSX file was wrong in v1, #462 +- corrected sign for spin-Y direction, PDOS, spin_moment, #486 - RealSpaceSI for right semi-infinite directions, #475 - tbtrans files now have a separate entry in the documentation ### Changed - changed DIIS solver to assume the matrix is symmetric (it is) - -### Changed - tbtncSileTBtrans and its derivates has changed, drastically. This will accommodate changes related to #477 and #478. Now `*_transmission` refers to energy resolved transmissions diff --git a/sisl/physics/electron.py b/sisl/physics/electron.py index 2e49703afd..bca7120bcf 100644 --- a/sisl/physics/electron.py +++ b/sisl/physics/electron.py @@ -229,16 +229,22 @@ def dot(v): # Initialize data PDOS = np.empty([4, state.shape[1] // 2, len(E)], dtype=dtype_complex_to_real(state.dtype)) + # Do spin-box calculations: + # PDOS[0] = total DOS (diagonal) + # PDOS[1] = x == < psi | \sigma_x S | psi > + # PDOS[2] = y == < psi | \sigma_y S | psi > + # PDOS[3] = z == < psi | \sigma_z S | psi > + d = distribution(E - eig[0]).reshape(1, -1) cs = conj(state[0]).reshape(-1, 2) v = S.dot(state[0].reshape(-1, 2)) D1 = (cs * v).real # uu,dd PDOS PDOS[0, :, :] = D1.sum(1).reshape(-1, 1) * d # total DOS PDOS[3, :, :] = (D1[:, 0] - D1[:, 1]).reshape(-1, 1) * d # z-dos - D1 = (cs[:, 1] * v[:, 0]).reshape(-1, 1) - D2 = (cs[:, 0] * v[:, 1]).reshape(-1, 1) + D1 = (cs[:, 1] * v[:, 0]).reshape(-1, 1) # d,u + D2 = (cs[:, 0] * v[:, 1]).reshape(-1, 1) # u,d PDOS[1, :, :] = (D1.real + D2.real) * d # x-dos - PDOS[2, :, :] = (D1.imag - D2.imag) * d # y-dos + PDOS[2, :, :] = (D2.imag - D1.imag) * d # y-dos for i in range(1, len(eig)): d = distribution(E - eig[i]).reshape(1, -1) cs = conj(state[i]).reshape(-1, 2) @@ -249,7 +255,7 @@ def dot(v): D1 = (cs[:, 1] * v[:, 0]).reshape(-1, 1) D2 = (cs[:, 0] * v[:, 1]).reshape(-1, 1) PDOS[1, :, :] += (D1.real + D2.real) * d - PDOS[2, :, :] += (D1.imag - D2.imag) * d + PDOS[2, :, :] += (D2.imag - D1.imag) * d else: PDOS = (conj(state[0]) * S.dot(state[0])).real.reshape(-1, 1) \ @@ -403,6 +409,8 @@ def dot(v): if S.shape[1] == state.shape[1]: S = S[::2, ::2] + # see PDOS for details related to the spin-box calculations + if project: s = np.empty([state.shape[0], state.shape[1] // 2, 3], dtype=dtype_complex_to_real(state.dtype)) @@ -414,7 +422,7 @@ def dot(v): D1 = cs[:, 1] * Sstate[:, 0] D2 = cs[:, 0] * Sstate[:, 1] s[i, :, 0] = D1.real + D2.real - s[i, :, 1] = D1.imag - D2.imag + s[i, :, 1] = D2.imag - D1.imag else: s = np.empty([state.shape[0], 3], dtype=dtype_complex_to_real(state.dtype)) @@ -427,9 +435,9 @@ def dot(v): cs = conj(state[i]).reshape(-1, 2) Sstate = S.dot(state[i].reshape(-1, 2)) D = cs.T @ Sstate - s[i, 2] = (D[0, 0] - D[1, 1]).real - s[i, 0] = (D[1, 0] + D[0, 1]).real - s[i, 1] = (D[1, 0] - D[0, 1]).imag + s[i, 2] = D[0, 0].real - D[1, 1].real + s[i, 0] = D[1, 0].real + D[0, 1].real + s[i, 1] = D[0, 1].imag - D[1, 0].imag return s diff --git a/sisl/physics/tests/test_hamiltonian.py b/sisl/physics/tests/test_hamiltonian.py index 0aa668ea99..5731f54a66 100644 --- a/sisl/physics/tests/test_hamiltonian.py +++ b/sisl/physics/tests/test_hamiltonian.py @@ -1116,6 +1116,76 @@ def test_pdos4(self, setup): assert PDOS.dtype.kind == 'f' assert np.allclose(PDOS.sum(0), DOS) + def test_pdos_nc(self): + geom = Geometry([0] * 3) + H = Hamiltonian(geom, spin="nc") + spin = H.spin + # this should be Hermitian + H[0, 0] = np.array([1, 2, 3, 4]) + E = [0] + def dist(E, *args): + return np.ones(len(E)) + + # just get a fictional PDOS + es = H.eigenstate() + PDOS = es.PDOS(E, dist)[..., 0] + SM = es.spin_moment() + SMp = es.spin_moment(project=True) + + # now check with spin stuff + pdos = es.inner().real + assert np.allclose(PDOS[0, 0], pdos.sum()) + + pdos = es.inner(matrix=spin.X).real + assert np.allclose(PDOS[1, 0], pdos.sum()) + assert np.allclose(SM[:, 0], pdos) + assert np.allclose(SMp[:, :, 0].sum(-1), pdos) + + pdos = es.inner(matrix=spin.Y).real + assert np.allclose(PDOS[2, 0], pdos.sum()) + assert np.allclose(SM[:, 1], pdos) + assert np.allclose(SMp[:, :, 1].sum(-1), pdos) + + pdos = es.inner(matrix=spin.Z).real + assert np.allclose(PDOS[3, 0], pdos.sum()) + assert np.allclose(SM[:, 2], pdos) + assert np.allclose(SMp[:, :, 2].sum(-1), pdos) + + def test_pdos_so(self): + geom = Geometry([0] * 3) + H = Hamiltonian(geom, spin="soc") + spin = H.spin + # this should be Hermitian + H[0, 0] = np.array([1, 2, 3, 4, 0, 0, 3, -4]) + E = [0] + def dist(E, *args): + return np.ones(len(E)) + + # just get a fictional PDOS + es = H.eigenstate() + PDOS = es.PDOS(E, dist)[..., 0] + SM = es.spin_moment() + SMp = es.spin_moment(project=True) + + # now check with spin stuff + pdos = es.inner().real + assert np.allclose(PDOS[0, 0], pdos.sum()) + + pdos = es.inner(matrix=spin.X).real + assert np.allclose(PDOS[1, 0], pdos.sum()) + assert np.allclose(SM[:, 0], pdos) + assert np.allclose(SMp[:, :, 0].sum(-1), pdos) + + pdos = es.inner(matrix=spin.Y).real + assert np.allclose(PDOS[2, 0], pdos.sum()) + assert np.allclose(SM[:, 1], pdos) + assert np.allclose(SMp[:, :, 1].sum(-1), pdos) + + pdos = es.inner(matrix=spin.Z).real + assert np.allclose(PDOS[3, 0], pdos.sum()) + assert np.allclose(SM[:, 2], pdos) + assert np.allclose(SMp[:, :, 2].sum(-1), pdos) + def test_coop_against_pdos_nonortho(self, setup): HS = setup.HS.copy() HS.construct([(0.1, 1.5), ((0., 1.), (1., 0.1))])