diff --git a/examples/socket/singlepoint/run.py b/examples/socket/singlepoint/spin-paired/run.py
similarity index 100%
rename from examples/socket/singlepoint/run.py
rename to examples/socket/singlepoint/spin-paired/run.py
diff --git a/examples/socket/singlepoint/spin-polarized/run.py b/examples/socket/singlepoint/spin-polarized/run.py
new file mode 100644
index 00000000..d31855da
--- /dev/null
+++ b/examples/socket/singlepoint/spin-polarized/run.py
@@ -0,0 +1,35 @@
+# Using SPARC
+from sparc.calculator import SPARC
+from ase.build import molecule
+import numpy as np
+
+water = molecule('H2O', vacuum=7)
+water.pbc = [False,False,False]
+water.set_initial_magnetic_moments([0.1, 0.1, 0.1])
+
+calc_params = {
+    "EXCHANGE_CORRELATION": "GGA_PBE",
+    "KPOINT_GRID": [1,1,1],
+    "MESH_SPACING": 0.35,
+    "TOL_SCF": 0.0001,
+    "MAXIT_SCF": 100,
+    "ELEC_TEMP_TYPE": "fermi-dirac",
+    "ELEC_TEMP": 116,
+    "PRINT_RESTART_FQ": 10,
+    "PRINT_ATOMS": 1,
+    "PRINT_FORCES": 1,
+    "SPIN_TYP": 1,
+}
+
+with SPARC(use_socket=True, **calc_params) as calc:
+    water.calc = calc
+    print("Initial magnetic moments before calculation call:\n", water.get_initial_magnetic_moments())
+    # energy = water.get_potential_energy()
+    # forces = water.get_forces()
+    net_magmom = water.get_magnetic_moment()
+    magmoms = water.get_magnetic_moments()
+# print('***********************************************************************************************')
+# print("Energy: {}\nMax Force: {}".format(energy, np.max(abs(forces)) ) )
+print('*'*100)
+print("Net magnetic moment: ", net_magmom)
+print("Atomic magnetic moments:\n", magmoms)
\ No newline at end of file
diff --git a/sparc/io.py b/sparc/io.py
index 8fccc625..eefc6611 100644
--- a/sparc/io.py
+++ b/sparc/io.py
@@ -534,6 +534,12 @@ def _extract_static_results(self, raw_results, index=":"):
             if "forces" in static_results:
                 partial_results["forces"] = static_results["forces"][self.resort]
 
+            if "atomic_magnetization" in static_results:
+                partial_results["magmoms"] = static_results["atomic_magnetization"][self.resort]
+
+            if "net_magnetization" in static_results:
+                partial_results["magmom"] = static_results["net_magnetization"]
+
             if "stress" in static_results:
                 partial_results["stress"] = static_results["stress"]
 
diff --git a/sparc/sparc_parsers/atoms.py b/sparc/sparc_parsers/atoms.py
index e3b7236f..a1a509db 100644
--- a/sparc/sparc_parsers/atoms.py
+++ b/sparc/sparc_parsers/atoms.py
@@ -104,7 +104,8 @@ def atoms_to_dict(
             block_dict["COORD"] = pos
         if write_spin:
             # TODO: should we process atoms with already calculated magmoms?
-            block_dict["SPIN"] = p_atoms.get_initial_magnetic_moments()
+            n_atom = len(p_atoms)
+            block_dict["SPIN"] = p_atoms.get_initial_magnetic_moments().reshape(n_atom,-1)
         if write_relax:
             relax_this_block = relax_mask[start:end]
             block_dict["RELAX"] = relax_this_block
diff --git a/sparc/sparc_parsers/ion.py b/sparc/sparc_parsers/ion.py
index 8bdb9859..280eefb3 100644
--- a/sparc/sparc_parsers/ion.py
+++ b/sparc/sparc_parsers/ion.py
@@ -124,6 +124,7 @@ def _write_ion(
             "PSEUDO_POT",
             "COORD_FRAC",
             "COORD",
+            "SPIN",
             "RELAX",
         ]:
             val = block.get(key, None)
@@ -139,7 +140,7 @@ def _write_ion(
             # TODO: make sure 1 line is accepted
             # TODO: write pads to vector lines
             if (val_string.count("\n") > 0) or (
-                key in ["COORD_FRAC", "COORD", "RELAX"]
+                key in ["COORD_FRAC", "COORD", "RELAX", "SPIN"]
             ):
                 output = f"{key}:\n{val_string}\n"
             else:
diff --git a/sparc/sparc_parsers/out.py b/sparc/sparc_parsers/out.py
index 16d189ed..b3b2ec34 100644
--- a/sparc/sparc_parsers/out.py
+++ b/sparc/sparc_parsers/out.py
@@ -222,6 +222,9 @@ def _read_scfs(contents):
             elif unit == "sec":
                 converted_value = raw_value * 1
                 converted_unit = "sec"
+            elif unit == "Bohr magneton":
+                converted_value = raw_value
+                converted_unit = "Bohr magneton"
             else:
                 warn(f"Conversion for unit {unit} unknown! Treat as unit")
                 converted_value = raw_value
diff --git a/sparc/sparc_parsers/static.py b/sparc/sparc_parsers/static.py
index 0058b8be..bba410d6 100644
--- a/sparc/sparc_parsers/static.py
+++ b/sparc/sparc_parsers/static.py
@@ -75,6 +75,10 @@ def _read_static_block(raw_block):
         name = "free energy"
     elif "Atomic forces" in header_name:
         name = "forces"
+    elif "Net magnetization" in header_name:
+        name = "net_magnetization"
+    elif "Atomic magnetization" in header_name:
+        name = "atomic_magnetization"
     elif "Stress (GPa)" in header_name:
         name = "stress"
     elif "Stress equiv." in header_name:
@@ -149,6 +153,10 @@ def _read_static_step(step):
             value = raw_value * Hartree
         elif name == "forces":
             value = raw_value * Hartree / Bohr
+        elif name == "atomic_magnetization":
+            value = raw_value
+        elif name == "net_magnetization":
+            value = raw_value
         elif name == "stress":
             # Stress is in eV/Ang^3, may need to convert to Virial later when cell is known
             # For low-dimension stress info, use stress_equiv