diff --git a/examples-py/hmc-pure-gauge-fa-sym.log b/examples-py/hmc-pure-gauge-fa-sym.log
index ce475a37a..14526bb78 100644
--- a/examples-py/hmc-pure-gauge-fa-sym.log
+++ b/examples-py/hmc-pure-gauge-fa-sym.log
@@ -1,29 +1,7 @@
-CHECK: ('job_tag', 'test-4nt8')
 CHECK: ('total_site', (4, 4, 4, 8))
-CHECK: ('load_config_params', {'twist_boundary_at_boundary': [0.0, 0.0, 0.0, -0.5]})
-CHECK: ('fermion_params', {0: {0: {'M5': 1.8, 'boundary_phases': [1.0, 1.0, 1.0, 1.0], 'b': 1.5, 'c': 0.5, 'Ls': 8, 'mass': 0.01}, 1: {'M5': 1.8, 'boundary_phases': [1.0, 1.0, 1.0, 1.0], 'b': 1.5, 'c': 0.5, 'Ls': 8, 'mass': 0.01}, 2: {'M5': 1.8, 'boundary_phases': [1.0, 1.0, 1.0, 1.0], 'b': 1.5, 'c': 0.5, 'Ls': 10, 'mass': 0.01}}, 1: {0: {'M5': 1.8, 'boundary_phases': [1.0, 1.0, 1.0, 1.0], 'b': 1.5, 'c': 0.5, 'Ls': 8, 'mass': 0.04}, 1: {'M5': 1.8, 'boundary_phases': [1.0, 1.0, 1.0, 1.0], 'b': 1.5, 'c': 0.5, 'Ls': 8, 'mass': 0.04}, 2: {'M5': 1.8, 'boundary_phases': [1.0, 1.0, 1.0, 1.0], 'b': 1.5, 'c': 0.5, 'Ls': 10, 'mass': 0.04}}, 2: {0: {'M5': 1.8, 'boundary_phases': [1.0, 1.0, 1.0, 1.0], 'b': 1.5, 'c': 0.5, 'Ls': 8, 'mass': 0.2}, 1: {'M5': 1.8, 'boundary_phases': [1.0, 1.0, 1.0, 1.0], 'b': 1.5, 'c': 0.5, 'Ls': 8, 'mass': 0.2}, 2: {'M5': 1.8, 'boundary_phases': [1.0, 1.0, 1.0, 1.0], 'b': 1.5, 'c': 0.5, 'Ls': 10, 'mass': 0.2}}})
-CHECK: ('lanc_params', {0: {0: {'fermion_params': {'M5': 1.8, 'boundary_phases': [1.0, 1.0, 1.0, 1.0], 'b': 1.5, 'c': 0.5, 'Ls': 8, 'mass': 0.01}, 'pit_params': {'eps': 0.01, 'maxiter': 500, 'real': True}, 'cheby_params': {'low': 0.1, 'high': 5.5, 'order': 50}, 'irl_params': {'Nstop': 20, 'Nk': 25, 'Nm': 30, 'resid': 1e-08, 'betastp': 0.0, 'maxiter': 20, 'Nminres': 1}}}})
-CHECK: ('clanc_params', {0: {0: {'block': [2, 2, 2, 2], 'nbasis': 20, 'cheby_params': {'low': 0.2, 'high': 5.5, 'order': 50}, 'irl_params': {'Nstop': 30, 'Nk': 35, 'Nm': 40, 'resid': 1e-08, 'betastp': 0.0, 'maxiter': 20, 'Nminres': 1}, 'smoother_params': {'eps': 1e-08, 'maxiter': 10}, 'save_params': {'nsingle': 10, 'mpi': [1, 1, 1, 4]}}}})
-CHECK: ('gf_ape_smear_coef', 0.5)
-CHECK: ('gf_ape_smear_step', 30)
-CHECK: ('prop_smear_coef', 0.9375)
-CHECK: ('prop_smear_step', 10)
-CHECK: ('trajs', [1000, 1100, 1200, 1300])
-CHECK: ('field-selection-fsel-rate', 0.0625)
-CHECK: ('field-selection-psel-rate', 0.03125)
-CHECK: ('field-selection-fsel-psrc-prop-norm-threshold', 0.001)
-CHECK: ('n_points_psel', 6)
-CHECK: ('n_exact_wsrc', 2)
-CHECK: ('n_per_tslice_smear', 2)
-CHECK: ('prob_acc_1_smear', 0.25)
-CHECK: ('prob_acc_2_smear', 0.0625)
-CHECK: ('prob_acc_1_psrc', 0.25)
-CHECK: ('prob_acc_2_psrc', 0.0625)
-CHECK: ('n_rand_u1_fsel', 4)
-CHECK: ('prob_acc_1_rand_u1', 0.25)
-CHECK: ('prob_acc_2_rand_u1', 0.0625)
-CHECK: ('meson_tensor_tsep', 1)
-CHECK: ('meson_jwjj_threshold', 0.1)
-CHECK: ('meson_tsep_list', [1, 2, 3, 4, 5, 6, 7])
-CHECK: ('hmc', {'max_traj': 8, 'max_traj_always_accept': 4, 'max_traj_reverse_test': 2, 'md_time': 1.0, 'n_step': 10, 'beta': 2.13, 'c1': -0.331, 'save_traj_interval': 4, 'is_saving_topo_info': True, 'fa': {'mass_type': 'random', 'complete_refresh_interval': 2}})
+CHECK: ('hmc', {'max_traj': 8, 'max_traj_always_accept': 4, 'max_traj_reverse_test': 2, 'md_time': 1.0, 'n_step': 10, 'beta': 2.13, 'c1': -0.331, 'save_traj_interval': 4, 'is_saving_topo_info': True})
+CHECK: ('total_site', (4, 4, 4, 8))
+CHECK: ('hmc', {'max_traj': 8, 'max_traj_always_accept': 4, 'max_traj_reverse_test': 2, 'md_time': 1.0, 'n_step': 10, 'beta': 2.13, 'c1': -0.331, 'save_traj_interval': 4, 'is_saving_topo_info': True, 'fa': {'complete_refresh_interval': 2, 'mass_type': 'random', 'interval_list': [1, 2]}})
+CHECK: ('total_site', (4, 4, 4, 8))
+CHECK: ('hmc', {'max_traj': 8, 'max_traj_always_accept': 4, 'max_traj_reverse_test': 2, 'md_time': 1.0, 'n_step': 10, 'beta': 2.13, 'c1': -0.331, 'save_traj_interval': 4, 'is_saving_topo_info': True, 'fa': {'complete_refresh_interval': 2, 'mass_type': 'grid-2', 'mass_list': [1.0, 1.25, 1.5, 4.0], 'interval_list': [1, 1, 1, 2]}})
 CHECK: finished successfully.
diff --git a/examples-py/hmc-pure-gauge-fa-sym.log.json b/examples-py/hmc-pure-gauge-fa-sym.log.json
index 69227fa6b..7e4f25823 100644
--- a/examples-py/hmc-pure-gauge-fa-sym.log.json
+++ b/examples-py/hmc-pure-gauge-fa-sym.log.json
@@ -1,34 +1,98 @@
 [
  [
   "run_hmc: 1 plaq",
-  0.7573486969674719
+  0.7450943920367527
  ],
  [
   "run_hmc: 2 plaq",
-  0.7318413424213871
+  0.6766937771184737
  ],
  [
   "run_hmc: 3 plaq",
-  0.6720538838704095
+  0.6436265119193414
  ],
  [
   "run_hmc: 4 plaq",
-  0.6697002531010241
+  0.6206781160838983
  ],
  [
   "run_hmc: 5 plaq",
-  0.6338969307784792
+  0.5973990387209166
  ],
  [
   "run_hmc: 6 plaq",
-  0.6319724682743478
+  0.5846407720593919
  ],
  [
   "run_hmc: 7 plaq",
-  0.6319724682743478
+  0.5867134719750315
  ],
  [
   "run_hmc: 8 plaq",
-  0.6236327145088388
+  0.5881418797220908
+ ],
+ [
+  "run_hmc: 1 plaq",
+  0.7574303159184317
+ ],
+ [
+  "run_hmc: 2 plaq",
+  0.7162682961967344
+ ],
+ [
+  "run_hmc: 3 plaq",
+  0.6688024944326814
+ ],
+ [
+  "run_hmc: 4 plaq",
+  0.6473990260130362
+ ],
+ [
+  "run_hmc: 5 plaq",
+  0.63015992725797
+ ],
+ [
+  "run_hmc: 6 plaq",
+  0.6183370383925257
+ ],
+ [
+  "run_hmc: 7 plaq",
+  0.6056014047510285
+ ],
+ [
+  "run_hmc: 8 plaq",
+  0.5978031451858511
+ ],
+ [
+  "run_hmc: 1 plaq",
+  0.7537772772647717
+ ],
+ [
+  "run_hmc: 2 plaq",
+  0.6936329797422313
+ ],
+ [
+  "run_hmc: 3 plaq",
+  0.6609974647924229
+ ],
+ [
+  "run_hmc: 4 plaq",
+  0.6395844999109757
+ ],
+ [
+  "run_hmc: 5 plaq",
+  0.6181505390074267
+ ],
+ [
+  "run_hmc: 6 plaq",
+  0.6103332490181326
+ ],
+ [
+  "run_hmc: 7 plaq",
+  0.6007069792836985
+ ],
+ [
+  "run_hmc: 8 plaq",
+  0.5964854465056361
  ]
 ]
\ No newline at end of file
diff --git a/examples-py/hmc-pure-gauge-fa-sym.py b/examples-py/hmc-pure-gauge-fa-sym.py
index 8a52a07be..a72289707 100755
--- a/examples-py/hmc-pure-gauge-fa-sym.py
+++ b/examples-py/hmc-pure-gauge-fa-sym.py
@@ -35,6 +35,7 @@ def gf_evolve(gf, gm, gm_dual, mf, mf_dual, dt):
 
 @q.timer
 def gm_evolve_fg(gm, gm_dual, gf_init, mf, mf_dual, ga, fg_dt, dt):
+    fname = q.get_fname()
     geo = gf_init.geo
     gf = GaugeField(geo)
     gf @= gf_init
@@ -95,7 +96,9 @@ def run_hmc_traj(
     gm_dual0 = GaugeMomentum(geo)
     gm0 @= gm
     gm_dual0 @= gm_dual
-    project_gauge_transform(gm0, gm_dual0, mf, mf_dual)
+    change_qnorm = project_gauge_transform(gm0, gm_dual0, mf, mf_dual)
+    if change_qnorm > 1e-8:
+        q.displayln_info(f"{fname}: project_gauge_transform: change_qnorm={change_qnorm}")
     delta_h = run_hmc_evolve(gm0, gm_dual0, gf0, mf, mf_dual, ga, rs, n_step, md_time)
     if is_reverse_test:
         gm_r = GaugeMomentum(geo)
@@ -137,36 +140,106 @@ def run_hmc_mass_mom_refresh(
         rs,
         mf, mf_dual,
         gm, gm_dual,
+        sel_list,
+        sel_dual_list,
         ):
     fname = q.get_fname()
     rs = rs.split(f"{traj}")
-    complete_refresh_interval = get_param(job_tag, "hmc", "fa", "complete_refresh_interval")
+    total_site = q.Coordinate(get_param(job_tag, "total_site"))
+    geo = q.Geometry(total_site)
+    complete_refresh_interval = get_param(job_tag, "hmc", "fa", "complete_refresh_interval", default=1)
     mass_type = get_param(job_tag, "hmc", "fa", "mass_type")
+    mass_list = get_param(job_tag, "hmc", "fa", "mass_list")
+    interval_list = get_param(job_tag, "hmc", "fa", "interval_list")
     is_refresh = traj % complete_refresh_interval == 0
     if is_force_refresh:
         if not is_refresh:
             q.displayln_info(f"{fname}: Force complete refresh of mass and momentum. {job_tag} {traj}")
         is_refresh = True
     if is_refresh:
+        # Completely refresh all mass and momentum (and selection)
         if mass_type is None:
             q.set_unit(mf)
             q.set_unit(mf_dual)
+            sel = np.zeros(mf[:].shape, dtype=bool)
+            sel[:] = True
+            sel_list[:] = [
+                    sel,
+                    ]
+            sel_dual_list[:] = [
+                    sel,
+                    ]
         elif mass_type == "random":
             mf.set_rand(rs.split("fa_mass"), 4.0, 1.0)
             mf_dual.set_rand(rs.split("fa_mass_dual"), 4.0, 1.0)
+            sel = mf[:] < 2.0
+            sel_dual = mf_dual[:] < 2.0
+            sel_list[:] = [
+                    sel,
+                    ~sel,
+                    ]
+            sel_dual_list[:] = [
+                    sel_dual,
+                    ~sel_dual,
+                    ]
+            assert isinstance(interval_list, list)
+            assert len(sel_list) == len(interval_list)
+        elif mass_type == "grid-2":
+            xg_arr = geo.xg_arr()
+            c_offset = rs.c_rand_gen(total_site).to_numpy()
+            xg_rel_arr = xg_arr - c_offset
+            count = np.zeros(xg_arr.shape, np.int32)
+            for m in range(4):
+                for n in range(4):
+                    if n == m:
+                        continue
+                    count[xg_rel_arr[:, n] % 2 == 0, m] += 1
+            assert np.all(count >= 0)
+            assert np.all(count <= 3)
+            sel_list[:] = [
+                count == 0,
+                count == 1,
+                count == 2,
+                count == 3,
+            ]
+            sel_dual_list[:] = sel_list
+            assert isinstance(interval_list, list)
+            assert len(sel_list) == len(interval_list)
+            assert isinstance(mass_list, list)
+            assert len(mass_list) == len(interval_list)
+            for idx, mass in enumerate(mass_list):
+                mf[sel_list[idx]] = mass
+                mf_dual[sel_dual_list[idx]] = mass
         else:
             raise Exception(f"{fname}: mass_type={mass_type}")
         gm.set_rand_fa(mf, rs.split("set_rand_gauge_momentum"))
         gm_dual.set_rand_fa(mf_dual, rs.split("set_rand_gauge_momentum_dual"))
     else:
-        sel = mf[:] < 2.0
-        sel_dual = mf_dual[:] < 2.0
-        gm1 = gm.copy()
-        gm_dual1 = gm_dual.copy()
+        # Only refresh some selection momentum (do not change mass and selection)
+        # IMPORTANT: Need to first add the gauge transform component back
+        gm1 = GaugeMomentum(geo)
+        gm_dual1 = GaugeMomentum(geo)
+        gm1.set_rand_fa(mf, rs.split("set_rand_gauge_momentum_g"))
+        gm_dual1.set_rand_fa(mf_dual, rs.split("set_rand_gauge_momentum_dual_g"))
+        gm2 = gm1.copy()
+        gm_dual2 = gm_dual1.copy()
+        project_gauge_transform(gm2, gm_dual2, mf, mf_dual)
+        gm2 -= gm1
+        gm_dual2 -= gm_dual1
+        gm -= gm2
+        gm_dual -= gm_dual2
+        # Then update the selected momentum
         gm1.set_rand_fa(mf, rs.split("set_rand_gauge_momentum"))
         gm_dual1.set_rand_fa(mf_dual, rs.split("set_rand_gauge_momentum_dual"))
-        gm[sel] = gm1[sel]
-        gm_dual[sel] = gm_dual1[sel]
+        assert len(sel_list) == len(interval_list)
+        assert len(sel_dual_list) == len(interval_list)
+        for sel, interval in zip(sel_list, interval_list):
+            if traj % interval == 0:
+                gm[sel] = gm1[sel]
+        for sel_dual, interval in zip(sel_dual_list, interval_list):
+            if traj % interval == 0:
+                gm_dual[sel_dual] = gm_dual1[sel_dual]
+    # Project out the gauge transform movement
     project_gauge_transform(gm, gm_dual, mf, mf_dual)
 
 @q.timer_verbose
@@ -221,6 +294,8 @@ def run_hmc(job_tag):
         traj = traj_load
         gf.load(get_load_path(f"{job_tag}/configs/ckpoint_lat.{traj}"))
     is_force_refresh = True
+    sel_list = []
+    sel_dual_list = []
     for traj in range(traj, max_traj):
         run_hmc_mass_mom_refresh(
                 job_tag, traj,
@@ -228,6 +303,8 @@ def run_hmc(job_tag):
                 rs.split("run_hmc_mass_mom_refresh"),
                 mf, mf_dual,
                 gm, gm_dual,
+                sel_list,
+                sel_dual_list,
                 )
         is_force_refresh = False
         traj += 1
@@ -260,7 +337,19 @@ def run_hmc(job_tag):
 
 # ----
 
-job_tag = "test-4nt8"
+job_tag = "test0-4nt8"
+set_param(job_tag, "total_site")((4, 4, 4, 8,))
+set_param(job_tag, "hmc", "max_traj")(8)
+set_param(job_tag, "hmc", "max_traj_always_accept")(4)
+set_param(job_tag, "hmc", "max_traj_reverse_test")(2)
+set_param(job_tag, "hmc", "md_time")(1.0)
+set_param(job_tag, "hmc", "n_step")(10)
+set_param(job_tag, "hmc", "beta")(2.13)
+set_param(job_tag, "hmc", "c1")(-0.331)
+set_param(job_tag, "hmc", "save_traj_interval")(4)
+set_param(job_tag, "hmc", "is_saving_topo_info")(True)
+
+job_tag = "test1-4nt8"
 set_param(job_tag, "total_site")((4, 4, 4, 8,))
 set_param(job_tag, "hmc", "max_traj")(8)
 set_param(job_tag, "hmc", "max_traj_always_accept")(4)
@@ -271,8 +360,25 @@ def run_hmc(job_tag):
 set_param(job_tag, "hmc", "c1")(-0.331)
 set_param(job_tag, "hmc", "save_traj_interval")(4)
 set_param(job_tag, "hmc", "is_saving_topo_info")(True)
+set_param(job_tag, "hmc", "fa", "complete_refresh_interval")(2)
 set_param(job_tag, "hmc", "fa", "mass_type")("random")
+set_param(job_tag, "hmc", "fa", "interval_list")([ 1, 2, ])
+
+job_tag = "test2-4nt8"
+set_param(job_tag, "total_site")((4, 4, 4, 8,))
+set_param(job_tag, "hmc", "max_traj")(8)
+set_param(job_tag, "hmc", "max_traj_always_accept")(4)
+set_param(job_tag, "hmc", "max_traj_reverse_test")(2)
+set_param(job_tag, "hmc", "md_time")(1.0)
+set_param(job_tag, "hmc", "n_step")(10)
+set_param(job_tag, "hmc", "beta")(2.13)
+set_param(job_tag, "hmc", "c1")(-0.331)
+set_param(job_tag, "hmc", "save_traj_interval")(4)
+set_param(job_tag, "hmc", "is_saving_topo_info")(True)
 set_param(job_tag, "hmc", "fa", "complete_refresh_interval")(2)
+set_param(job_tag, "hmc", "fa", "mass_type")("grid-2")
+set_param(job_tag, "hmc", "fa", "mass_list")([ 1.0, 1.25, 1.5, 4.0, ])
+set_param(job_tag, "hmc", "fa", "interval_list")([ 1, 1, 1, 2, ])
 
 job_tag = "32I_b2p8_fa_sym"
 set_param(job_tag, "total_site")((32, 32, 32, 64,))
@@ -339,6 +445,38 @@ def run_hmc(job_tag):
 set_param(job_tag, "hmc", "save_traj_interval")(2)
 set_param(job_tag, "hmc", "is_saving_topo_info")(True)
 
+job_tag = "32I_b2p8_fa_sym_v1"
+set_param(job_tag, "total_site")((32, 32, 32, 64,))
+set_param(job_tag, "a_inv_gev")(2.646) # 2003 lattice spacing 0309017.pdf
+set_param(job_tag, "hmc", "max_traj")(20000)
+set_param(job_tag, "hmc", "max_traj_always_accept")(100)
+set_param(job_tag, "hmc", "max_traj_reverse_test")(2)
+set_param(job_tag, "hmc", "md_time")(1.0)
+set_param(job_tag, "hmc", "n_step")(32)
+set_param(job_tag, "hmc", "beta")(2.80)
+set_param(job_tag, "hmc", "c1")(-0.331)
+set_param(job_tag, "hmc", "save_traj_interval")(10)
+set_param(job_tag, "hmc", "is_saving_topo_info")(True)
+set_param(job_tag, "hmc", "fa", "mass_type")("grid-2")
+set_param(job_tag, "hmc", "fa", "mass_list")([ 1.0, 1.25, 1.5, 4.0, ])
+set_param(job_tag, "hmc", "fa", "interval_list")([ 1, 1, 1, 2, ])
+
+job_tag = "32I_b2p8_fa_sym_v2"
+set_param(job_tag, "total_site")((32, 32, 32, 64,))
+set_param(job_tag, "a_inv_gev")(2.646) # 2003 lattice spacing 0309017.pdf
+set_param(job_tag, "hmc", "max_traj")(20000)
+set_param(job_tag, "hmc", "max_traj_always_accept")(100)
+set_param(job_tag, "hmc", "max_traj_reverse_test")(2)
+set_param(job_tag, "hmc", "md_time")(2.0)
+set_param(job_tag, "hmc", "n_step")(32)
+set_param(job_tag, "hmc", "beta")(2.80)
+set_param(job_tag, "hmc", "c1")(-0.331)
+set_param(job_tag, "hmc", "save_traj_interval")(10)
+set_param(job_tag, "hmc", "is_saving_topo_info")(True)
+set_param(job_tag, "hmc", "fa", "mass_type")("grid-2")
+set_param(job_tag, "hmc", "fa", "mass_list")([ 1.0, 1.25, 1.5, 4.0, ])
+set_param(job_tag, "hmc", "fa", "interval_list")([ 1, 1, 1, 2, ])
+
 # ----
 
 size_node_list = [
@@ -364,7 +502,9 @@ def run_hmc(job_tag):
     #######################################################
 
     job_tags_default = [
-            "test-4nt8",
+            "test0-4nt8",
+            "test1-4nt8",
+            "test2-4nt8",
             ]
 
     if job_tags == [ "", ]:
diff --git a/qlat/qlat/everything.pxd b/qlat/qlat/everything.pxd
index b9043944e..00b8edc43 100644
--- a/qlat/qlat/everything.pxd
+++ b/qlat/qlat/everything.pxd
@@ -338,8 +338,7 @@ cdef extern from "qlat/hmc.h" namespace "qlat":
     void gf_evolve_dual(GaugeField& gf, const GaugeMomentum& gm_dual, const Field[RealD]& mf_dual, const RealD step_size) except +
     void set_gm_force(GaugeMomentum& gm_force, const GaugeField& gf, const GaugeAction& ga) except +
     void set_gm_force_dual(GaugeMomentum& gm_force_dual, const GaugeField& gf, const GaugeMomentum& gm_force) except +
-    void project_gauge_transform(GaugeMomentum& gm, GaugeMomentum& gm_dual, const Field[RealD]& mf, const Field[RealD]& mf_dual) except +
-
+    RealD project_gauge_transform(GaugeMomentum& gm, GaugeMomentum& gm_dual, const Field[RealD]& mf, const Field[RealD]& mf_dual) except +
 
 cdef extern from "qlat/vector_utils/utils_smear_vecs.h" namespace "qlat":
 
diff --git a/qlat/qlat/include/qlat/hmc.h b/qlat/qlat/include/qlat/hmc.h
index bb3b1ffdc..4aa1ed790 100644
--- a/qlat/qlat/include/qlat/hmc.h
+++ b/qlat/qlat/include/qlat/hmc.h
@@ -60,9 +60,9 @@ void set_gm_force(GaugeMomentum& gm_force, const GaugeField& gf,
 void set_gm_force_dual(GaugeMomentum& gm_force_dual, const GaugeField& gf,
                        const GaugeMomentum& gm_force);
 
-void project_gauge_transform(GaugeMomentum& gm, GaugeMomentum& gm_dual,
-                             const Field<RealD>& mf,
-                             const Field<RealD>& mf_dual);
+RealD project_gauge_transform(GaugeMomentum& gm, GaugeMomentum& gm_dual,
+                              const Field<RealD>& mf,
+                              const Field<RealD>& mf_dual);
 
 // -------------------
 
diff --git a/qlat/qlat/lib/hmc.cpp b/qlat/qlat/lib/hmc.cpp
index 128d43a3a..8676b32ce 100644
--- a/qlat/qlat/lib/hmc.cpp
+++ b/qlat/qlat/lib/hmc.cpp
@@ -369,9 +369,11 @@ void set_gm_force_dual(GaugeMomentum& gm_force_dual, const GaugeField& gf,
   });
 }
 
-void project_gauge_transform(GaugeMomentum& gm, GaugeMomentum& gm_dual,
-                             const Field<RealD>& mf,
-                             const Field<RealD>& mf_dual)
+RealD project_gauge_transform(GaugeMomentum& gm, GaugeMomentum& gm_dual,
+                              const Field<RealD>& mf,
+                              const Field<RealD>& mf_dual)
+// return qnorm of the change / geo.total_volume()
+// (up to an overall coefficient)
 {
   TIMER("project_gauge_transform");
   const Geometry& geo = gm.geo();
@@ -401,7 +403,6 @@ void project_gauge_transform(GaugeMomentum& gm, GaugeMomentum& gm_dual,
     const Coordinate xl = geo.coordinate_from_index(index);
     const Vector<ColorMatrix> v1 = gm.get_elems_const(xl);
     const Vector<RealD> vm1 = mf.get_elems_const(xl);
-    RealD inv_m_sum_x = 0.0;
     ColorMatrix v_sum;
     set_zero(v_sum);
     for (int m = 0; m < 4; ++m) {
@@ -410,12 +411,10 @@ void project_gauge_transform(GaugeMomentum& gm, GaugeMomentum& gm_dual,
       const RealD m2m = mf2_ext.get_elem(xl_m, m);
       const ColorMatrix& c1 = v1[m];
       const ColorMatrix& c2m = gm2_ext.get_elem(xl_m, m);
-      inv_m_sum_x += 1.0 / m1 + 1.0 / m2m;
-      v_sum += c1 * (1.0 / m1) + c2m * (1.0 / m2m);
+      v_sum += c1 * (1.0 / std::sqrt(m1)) + c2m * (1.0 / std::sqrt(m2m));
     }
-    const RealD alpha = 8.0 / inv_m_sum_x;
     ColorMatrix& ag = gm_ag_ext.get_elem(xl);
-    ag = (alpha / 8.0) * v_sum;
+    ag = (1.0 / 8.0) * v_sum;
   });
   refresh_expanded_1(gm_ag_ext);
   qacc_for(index, geo.local_volume(), {
@@ -423,13 +422,16 @@ void project_gauge_transform(GaugeMomentum& gm, GaugeMomentum& gm_dual,
     const Coordinate xl = geo.coordinate_from_index(index);
     Vector<ColorMatrix> v1 = gm.get_elems(index);
     Vector<ColorMatrix> v2 = gm_dual.get_elems(index);
+    const Vector<RealD> vm1 = mf.get_elems_const(xl);
+    const Vector<RealD> vm2 = mf_dual.get_elems_const(xl);
     const ColorMatrix& ag1 = gm_ag_ext.get_elem(xl);
     for (int m = 0; m < 4; ++m) {
       const Coordinate xl_p = coordinate_shifts(xl, m);
-      v1[m] -= ag1;
-      v2[m] -= gm_ag_ext.get_elem(xl_p);
+      v1[m] -= std::sqrt(vm1[m]) * ag1;
+      v2[m] -= std::sqrt(vm2[m]) * gm_ag_ext.get_elem(xl_p);
     }
   });
+  return qnorm(gm_ag_ext) / geo.total_volume();
 }
 
 }  // namespace qlat