Fix issue 486 - in correct sign of the y component of the spin moments

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

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
 

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))])