diff --git a/desc/input_reader.py b/desc/input_reader.py
index 2316c8bd5b..e59a17eba0 100644
--- a/desc/input_reader.py
+++ b/desc/input_reader.py
@@ -1335,7 +1335,7 @@ def parse_vmec_inputs(vmec_fname, threshold=0):  # noqa: C901 - FIXME: simplify
                         )
                         inputs["axis"][-1, :] = np.array([n, -numbers[k], 0.0])
             match = re.search(
-                r"ZAXIS_CC\s*=(\s*" + num_form + r"\s*,?)*", command, re.IGNORECASE
+                r"ZAXIS(_CS)?\s*=(\s*" + num_form + r"\s*,?)*", command, re.IGNORECASE
             )
             if match:
                 numbers = [
@@ -1343,18 +1343,19 @@ def parse_vmec_inputs(vmec_fname, threshold=0):  # noqa: C901 - FIXME: simplify
                     for x in re.findall(num_form, match.group(0))
                     if re.search(r"\d", x)
                 ]
-                for k in range(len(numbers)):
-                    n = k
+                # ignore the n=0 since it should always be zero for sin terms
+                for k in range(1, len(numbers)):
+                    n = -k
                     idx = np.where(inputs["axis"][:, 0] == n)[0]
                     if np.size(idx):
-                        inputs["axis"][idx[0], 2] = numbers[k]
+                        inputs["axis"][idx[0], 2] = -numbers[k]
                     else:
                         inputs["axis"] = np.pad(
                             inputs["axis"], ((0, 1), (0, 0)), mode="constant"
                         )
-                        inputs["axis"][-1, :] = np.array([n, 0.0, numbers[k]])
+                        inputs["axis"][-1, :] = np.array([n, 0.0, -numbers[k]])
             match = re.search(
-                r"ZAXIS_(CS)?\s*=(\s*" + num_form + r"\s*,?)*", command, re.IGNORECASE
+                r"ZAXIS_CC\s*=(\s*" + num_form + r"\s*,?)*", command, re.IGNORECASE
             )
             if match:
                 numbers = [
@@ -1362,17 +1363,16 @@ def parse_vmec_inputs(vmec_fname, threshold=0):  # noqa: C901 - FIXME: simplify
                     for x in re.findall(num_form, match.group(0))
                     if re.search(r"\d", x)
                 ]
-                # ignore the n=0 since it should always be zero for sin terms
-                for k in range(1, len(numbers)):
-                    n = -k
+                for k in range(len(numbers)):
+                    n = k
                     idx = np.where(inputs["axis"][:, 0] == n)[0]
                     if np.size(idx):
-                        inputs["axis"][idx[0], 2] = -numbers[k]
+                        inputs["axis"][idx[0], 2] = numbers[k]
                     else:
                         inputs["axis"] = np.pad(
                             inputs["axis"], ((0, 1), (0, 0)), mode="constant"
                         )
-                        inputs["axis"][-1, :] = np.array([n, 0.0, -numbers[k]])
+                        inputs["axis"][-1, :] = np.array([n, 0.0, numbers[k]])
 
             # boundary shape
             # RBS*sin(m*t-n*p) = RBS*sin(m*t)*cos(n*p) - RBS*cos(m*t)*sin(n*p)