diff --git a/CHANGELOG.md b/CHANGELOG.md
index bb801a909..b2c1f9ee9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -87,6 +87,7 @@ we hit release version 1.0.0.
 
 
 ### Fixed
+- fixed Mulliken calculations for polarized calculations due to missing copy, #611
 - fixed single argument `ret_isc=True` of `close`, #604 and #605
 - tiling Grid now only possible for commensurate grids (grid.lattice % grid.geometry.lattice)
 - rare cases for non-Gamma calculations with actual Gamma matrices resulted
diff --git a/src/sisl/physics/densitymatrix.py b/src/sisl/physics/densitymatrix.py
index 923502602..79b81e0d5 100644
--- a/src/sisl/physics/densitymatrix.py
+++ b/src/sisl/physics/densitymatrix.py
@@ -326,7 +326,8 @@ def _convert(M):
                 M[3] = 0.5 * (M[3] - M[7]) # sign change again below
                 M = M[:4]
             elif M.shape[0] == 2:
-                tmp = M[1]
+                # necessary to not overwrite data
+                tmp = M[1].copy()
                 M[1] = M[0] - M[1]
                 M[0] += tmp
             elif M.shape[0] == 1:
diff --git a/src/sisl/physics/tests/test_density_matrix.py b/src/sisl/physics/tests/test_density_matrix.py
index c8363a31c..5c54e9ab2 100644
--- a/src/sisl/physics/tests/test_density_matrix.py
+++ b/src/sisl/physics/tests/test_density_matrix.py
@@ -22,7 +22,7 @@
 def setup():
     class t():
         def __init__(self):
-            bond = 1.42
+            self.bond = bond = 1.42
             sq3h = 3.**.5 * 0.5
             self.lattice = Lattice(np.array([[1.5, sq3h, 0.],
                                              [1.5, -sq3h, 0.],
@@ -127,6 +127,30 @@ def test_mulliken_values_non_orthogonal(self, setup):
         assert mulliken[0] == pytest.approx(4)
         assert mulliken.sum() == pytest.approx(4)
 
+    def test_mulliken_polarized(self):
+        bond = 1.42
+        sq3h = 3.**.5 * 0.5
+        lattice = Lattice(np.array([[1.5, sq3h, 0.],
+                                 [1.5, -sq3h, 0.],
+                                 [0., 0., 10.]], np.float64) * bond, nsc=[3, 3, 1])
+
+        orb = AtomicOrbital('px', R=bond * 1.001)
+        C = Atom(6, orb)
+        g = Geometry(np.array([[0., 0., 0.],
+                               [1., 0., 0.]], np.float64) * bond,
+                     atoms=C, lattice=lattice)
+        D = DensityMatrix(g, spin=Spin('P'))
+        # 1 charge onsite for each spin-up
+        # 0.5 charge onsite for each spin-down
+        D.construct([[0.1, bond + 0.01], [(1., 0.5), (0.1, 0.1)]])
+
+        m = D.mulliken("orbital")
+        assert m[0].sum() == pytest.approx(3)
+        assert m[1].sum() == pytest.approx(1)
+        m = D.mulliken("atom")
+        assert m[0].sum() == pytest.approx(3)
+        assert m[1].sum() == pytest.approx(1)
+
     def test_rho1(self, setup):
         D = setup.D.copy()
         D.construct(setup.func)
@@ -134,7 +158,7 @@ def test_rho1(self, setup):
         D.density(grid)
 
     @pytest.mark.filterwarnings("ignore", message="*is NOT Hermitian for on-site")
-    def test_rho2(self, setup):
+    def test_rho2(self):
         bond = 1.42
         sq3h = 3.**.5 * 0.5
         lattice = Lattice(np.array([[1.5, sq3h, 0.],
@@ -178,7 +202,7 @@ def test_rho2(self, setup):
         D.density(grid, Spin.Z)
 
     @pytest.mark.filterwarnings("ignore", message="*is NOT Hermitian for on-site")
-    def test_orbital_momentum(self, setup):
+    def test_orbital_momentum(self):
         bond = 1.42
         sq3h = 3.**.5 * 0.5
         lattice = Lattice(np.array([[1.5, sq3h, 0.],
@@ -195,7 +219,7 @@ def test_orbital_momentum(self, setup):
         D.orbital_momentum("atom")
         D.orbital_momentum("orbital")
 
-    def test_spin_align_pol(self, setup):
+    def test_spin_align_pol(self):
         bond = 1.42
         sq3h = 3.**.5 * 0.5
         lattice = Lattice(np.array([[1.5, sq3h, 0.],
@@ -220,7 +244,7 @@ def test_spin_align_pol(self, setup):
         assert not np.allclose(D_mull[1], d_mull[3])
         assert np.allclose(D_mull[0], d_mull[0])
 
-    def test_spin_align_nc(self, setup):
+    def test_spin_align_nc(self):
         bond = 1.42
         sq3h = 3.**.5 * 0.5
         lattice = Lattice(np.array([[1.5, sq3h, 0.],
@@ -242,7 +266,7 @@ def test_spin_align_nc(self, setup):
         assert np.allclose(D_mull[0], d_mull[0])
 
     @pytest.mark.filterwarnings("ignore", message="*is NOT Hermitian for on-site")
-    def test_spin_align_so(self, setup):
+    def test_spin_align_so(self):
         bond = 1.42
         sq3h = 3.**.5 * 0.5
         lattice = Lattice(np.array([[1.5, sq3h, 0.],
@@ -263,7 +287,7 @@ def test_spin_align_so(self, setup):
         assert not np.allclose(D_mull, d_mull)
         assert np.allclose(D_mull[0], d_mull[0])
 
-    def test_spin_rotate_pol(self, setup):
+    def test_spin_rotate_pol(self):
         bond = 1.42
         sq3h = 3.**.5 * 0.5
         lattice = Lattice(np.array([[1.5, sq3h, 0.],
@@ -288,7 +312,7 @@ def test_spin_rotate_pol(self, setup):
         assert not np.allclose(D_mull[1], d_mull[3])
         assert np.allclose(D_mull[0], d_mull[0])
 
-    def test_spin_rotate_nc(self, setup):
+    def test_spin_rotate_nc(self):
         bond = 1.42
         sq3h = 3.**.5 * 0.5
         lattice = Lattice(np.array([[1.5, sq3h, 0.],
@@ -312,7 +336,7 @@ def test_spin_rotate_nc(self, setup):
         assert np.allclose(D_mull[0], d_mull[0])
 
     @pytest.mark.filterwarnings("ignore", message="*is NOT Hermitian for on-site")
-    def test_spin_rotate_so(self, setup):
+    def test_spin_rotate_so(self):
         bond = 1.42
         sq3h = 3.**.5 * 0.5
         lattice = Lattice(np.array([[1.5, sq3h, 0.],
@@ -345,7 +369,7 @@ def test_rho_smaller_grid1(self, setup):
         grid = Grid(0.2, geometry=setup.D.geometry.copy(), lattice=lattice)
         D.density(grid)
 
-    def test_rho_fail_p(self, setup):
+    def test_rho_fail_p(self):
         bond = 1.42
         sq3h = 3.**.5 * 0.5
         lattice = Lattice(np.array([[1.5, sq3h, 0.],
@@ -367,7 +391,7 @@ def test_rho_fail_p(self, setup):
         with pytest.raises(ValueError):
             D.density(grid, [1., -1, 0.])
 
-    def test_rho_fail_nc(self, setup):
+    def test_rho_fail_nc(self):
         bond = 1.42
         sq3h = 3.**.5 * 0.5
         lattice = Lattice(np.array([[1.5, sq3h, 0.],