diff --git a/examples/3_Advanced/coil_force_optimization/analysis_tools.py b/examples/3_Advanced/coil_force_optimization/analysis_tools.py
index 989d70342..6cdb205e6 100644
--- a/examples/3_Advanced/coil_force_optimization/analysis_tools.py
+++ b/examples/3_Advanced/coil_force_optimization/analysis_tools.py
@@ -70,9 +70,9 @@ def get_dfs(INPUT_DIR='./output/QA/with-force-penalty/1/optimizations/', OUTPUT_
         f"and max_arclength_variance < 1e-2"
         f"and coil_surface_distance < 0.375"
         f"and coil_coil_distance < 0.15"
-        f"and max_length > 4.0"
-        f"and normalized_BdotN < {2e-2}"
-        f"and max_max_force<50000"
+        f"and max_length > 3.5"
+        f"and normalized_BdotN < {4e-2}"
+        f"and max_max_force<60000"
     )
 
     ### STEP 3: Generate Pareto front and export UUIDs as .txt
diff --git a/examples/3_Advanced/coil_force_optimization/coil_force_pareto_scans.py b/examples/3_Advanced/coil_force_optimization/coil_force_pareto_scans.py
index 5a2553644..710579ebb 100644
--- a/examples/3_Advanced/coil_force_optimization/coil_force_pareto_scans.py
+++ b/examples/3_Advanced/coil_force_optimization/coil_force_pareto_scans.py
@@ -8,16 +8,16 @@
 # import imageio
 import matplotlib.pyplot as plt
 
-INPUT_DIR = "./output/QA/TVE/"
+INPUT_DIR = "./output/QA/TVE_new/"
 df, df_filtered, df_pareto = get_dfs(INPUT_DIR=INPUT_DIR)
 success_plt(df, df_filtered).show()
 df, df_filtered, df_pareto = get_dfs(INPUT_DIR=INPUT_DIR)
 
-df_filtered = df
-df_pareto = df
+#df_filtered = df
+#df_pareto = df
 y_axes = ["tve"]  #, "mean_RMS_force", "max_max_torque", "mean_RMS_torque", "net_forces", "net_torques"]
 labels = ["Total Vacuum Magnetic Energy [J]"]  #, "mean force [N/m]", "max torque [N]", "mean torque [N]", "net force [N]", "net torque [N-m]"]
-y_lims = [(10000, 35000)]  #, (5500, 10000), (5000, 30000), (2000, 6200), (4000, 13000), (5e2, 7000)]
+y_lims = [(7500, 30000)]  #, (5500, 10000), (5000, 30000), (2000, 6200), (4000, 13000), (5e2, 7000)]
 #keys = ["coil_coil_distance"]
 #keys = ["max_length"]
 #keys = ["max_max_κ"]
@@ -102,8 +102,8 @@
         )
         plt.xlabel(r'$\langle|\mathbf{B}\cdot\mathbf{n}|\rangle/\langle B \rangle$ [unitless]')
         plt.ylabel(label)
-        plt.xlim(0.7 * min(df_filtered["normalized_BdotN"]), max(df_filtered["normalized_BdotN"]))
-        #plt.ylim(y_lim)
+        plt.xlim(0.5 * min(df_filtered["normalized_BdotN"]), max(df_filtered["normalized_BdotN"]))
+        plt.ylim(y_lim)
         plt.xscale("log")
         #if y_lim[0] < 1:
         #    plt.yscale("log")
diff --git a/examples/3_Advanced/coil_force_optimization/cold_starts.sh b/examples/3_Advanced/coil_force_optimization/cold_starts.sh
index 28e2fb827..954107b47 100644
--- a/examples/3_Advanced/coil_force_optimization/cold_starts.sh
+++ b/examples/3_Advanced/coil_force_optimization/cold_starts.sh
@@ -1,12 +1,12 @@
 #!/bin/bash
-#SBATCH --time=600
+#SBATCH --time=500
 #SBATCH --mail-type=END
 #SBATCH --mail-user=aak572@nyu.edu
 #SBATCH --nodes=1
 #SBATCH --mem=48000
 #SBATCH --ntasks-per-node=32
 #SBATCH --cpus-per-task=1
-#SBATCH --array=0-30
+#SBATCH --array=0-20
 #SBATCH --partition=cs
 
 export OMP_NUM_THREADS=1
diff --git a/examples/3_Advanced/coil_force_optimization/initiation.py b/examples/3_Advanced/coil_force_optimization/initiation.py
index 15b574144..4884aca34 100644
--- a/examples/3_Advanced/coil_force_optimization/initiation.py
+++ b/examples/3_Advanced/coil_force_optimization/initiation.py
@@ -1,5 +1,5 @@
 from optimization_tools import *
 
-OUTPUT_DIR="./output/QA/TVE/"
+OUTPUT_DIR="./output/QA/TVE_fix/"
 INPUT_FILE="../../../tests/test_files/input.LandremanPaul2021_QA"
-initial_optimizations(N=20, OUTPUT_DIR=OUTPUT_DIR, INPUT_FILE=INPUT_FILE, with_force=True, FORCE_OBJ=TVE, MAXITER=100, debug=True)
+initial_optimizations(OUTPUT_DIR=OUTPUT_DIR, INPUT_FILE=INPUT_FILE, with_force=True, FORCE_OBJ=TVE)
diff --git a/src/simsopt/field/force.py b/src/simsopt/field/force.py
index 15c76c84e..e9c2c7b5d 100644
--- a/src/simsopt/field/force.py
+++ b/src/simsopt/field/force.py
@@ -1926,7 +1926,7 @@ def tve_pure(gamma, gammadash, gammas2, gammadashs2, current, currents2,
         gamma, gammadash, gammas2, gammadashs2, 
         quadpoints, quadpoints2, a, b, downsample, cross_section
     )
-    U = 0.5 * (jnp.sum(jnp.abs(Ii_Ij) * Lij[1:]) + Lij[0] * current ** 2)
+    U = 0.5 * (jnp.sum(Ii_Ij * Lij[1:]) + Lij[0] * current ** 2)
     return U
 
 class TVE(Optimizable):