diff --git a/doc/src/fix_electrode.rst b/doc/src/fix_electrode.rst
index 3d543f08d2..3759ff23ca 100644
--- a/doc/src/fix_electrode.rst
+++ b/doc/src/fix_electrode.rst
@@ -45,7 +45,7 @@ Syntax
                 rng_v = integer used to initialize random number generator
 
 * zero or more keyword/value pairs may be appended
-* keyword = *algo* or *symm* or *couple* or *etypes* or *ffield* or *write_mat* or *write_inv* or *read_mat* or *read_inv*
+* keyword = *algo* or *symm* or *couple* or *etypes* or *ffield* or *write_mat* or *write_inv* or *read_mat* or *read_inv* or *qtotal* or *eta*
 
 .. parsed-literal::
 
@@ -68,6 +68,10 @@ Syntax
         filename = file from which to read elastance matrix
     *read_inv* value = filename
         filename = file from which to read inverted matrix
+    *qtotal* value = number or *v_* equal-style variable
+        add overall potential so that all electrode charges add up to *qtotal*
+    *eta* value = d_propname
+        d_propname = a custom double vector defined via fix property/atom
 
 Examples
 """"""""
@@ -249,6 +253,24 @@ be enabled if any electrode particle has the same type as any
 electrolyte particle (which would be unusual in a typical simulation)
 and the fix will issue an error in that case.
 
+.. versionchanged:: qtotal
+
+The keyword *qtotal* causes *fix electrode/conp* and *fix electrode/thermo*
+to add an overall potential to all electrodes so that the total charge on
+the electrodes is a specified amount (which may be an equal-style variable).
+For example, if a user wanted to simulate a solution of excess cations
+such that the total electrolyte charge is +2, setting *qtotal -2* would cause
+the total electrode charge to be -2, so that the simulation box remains overall
+electroneutral. Since *fix electrode/conq* constrains the total charges of
+individual electrodes, and since *symm on* constrains the total charge of all
+electrodes to be zero, either option is incompatible with the *qtotal* keyword
+(even if *qtotal* is set to zero).
+
+The keyword *eta* takes the name of a custom double vector defined via fix
+property/atom.  The values will be used instead of the standard eta value.  The
+property/atom fix must be for vector of double values and use the *ghost on*
+option.
+
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
diff --git a/examples/PACKAGES/electrode/madelung/eval.py b/examples/PACKAGES/electrode/madelung/eval.py
index 2f5a355d9b..feda0e384e 100644
--- a/examples/PACKAGES/electrode/madelung/eval.py
+++ b/examples/PACKAGES/electrode/madelung/eval.py
@@ -1,7 +1,7 @@
 #!/usr/env/python3
 
-import sys
 import os.path as op
+import sys
 
 
 def rel_error(out, ref):
@@ -49,5 +49,5 @@ for label, ref, out in out_lines:
     error = rel_error(out, ref)
     lines.append(f"{label}: {out:.5f}, {error:.5f}\n")
 
-with open("madelung.txt", 'a') as f:
+with open("madelung.txt", "a") as f:
     f.writelines(lines)
diff --git a/examples/PACKAGES/electrode/madelung/in.eta b/examples/PACKAGES/electrode/madelung/in.eta
new file mode 100644
index 0000000000..3a45bb1bf5
--- /dev/null
+++ b/examples/PACKAGES/electrode/madelung/in.eta
@@ -0,0 +1,14 @@
+boundary p p f
+kspace_style ewald/electrode 1.0e-8
+kspace_modify slab 8.0 # ew3dc
+
+include "settings.mod" # styles, computes, groups and fixes
+
+thermo_style custom step pe c_qbot c_qtop
+fix feta all property/atom d_eta ghost on
+set group bot d_eta 2.0
+set group top d_eta 2.0
+fix conp bot electrode/conp 0 2 couple top 1 symm on eta d_eta write_inv inv.csv write_vec vec.csv
+
+run 0
+
diff --git a/examples/PACKAGES/electrode/madelung/in.eta_cg b/examples/PACKAGES/electrode/madelung/in.eta_cg
new file mode 100644
index 0000000000..5ac8cddf17
--- /dev/null
+++ b/examples/PACKAGES/electrode/madelung/in.eta_cg
@@ -0,0 +1,14 @@
+boundary p p f
+kspace_style ewald/electrode 1.0e-8
+kspace_modify slab 8.0 # ew3dc
+
+include "settings.mod" # styles, computes, groups and fixes
+
+thermo_style custom step pe c_qbot c_qtop
+fix feta all property/atom d_eta ghost on
+set group bot d_eta 0.5
+set group top d_eta 3.0
+fix conp bot electrode/conp 0 2 couple top 1 symm on eta d_eta algo cg 1e-6
+
+run 0
+
diff --git a/examples/PACKAGES/electrode/madelung/in.eta_mix b/examples/PACKAGES/electrode/madelung/in.eta_mix
new file mode 100644
index 0000000000..d00e008fa4
--- /dev/null
+++ b/examples/PACKAGES/electrode/madelung/in.eta_mix
@@ -0,0 +1,14 @@
+boundary p p f
+kspace_style ewald/electrode 1.0e-8
+kspace_modify slab 8.0 # ew3dc
+
+include "settings.mod" # styles, computes, groups and fixes
+
+thermo_style custom step pe c_qbot c_qtop
+fix feta all property/atom d_eta ghost on
+set group bot d_eta 0.5
+set group top d_eta 3.0
+fix conp bot electrode/conp 0 2 couple top 1 symm on eta d_eta write_inv inv.csv write_vec vec.csv
+
+run 0
+
diff --git a/examples/PACKAGES/electrode/madelung/log.19Feb2024.eta.g++.1 b/examples/PACKAGES/electrode/madelung/log.19Feb2024.eta.g++.1
new file mode 100644
index 0000000000..daf0563799
--- /dev/null
+++ b/examples/PACKAGES/electrode/madelung/log.19Feb2024.eta.g++.1
@@ -0,0 +1,138 @@
+LAMMPS (21 Nov 2023 - Development - patch_21Nov2023-668-g5b6c0c6b56)
+  using 1 OpenMP thread(s) per MPI task
+boundary p p f
+kspace_style ewald/electrode 1.0e-8
+kspace_modify slab 8.0 # ew3dc
+
+include "settings.mod" # styles, computes, groups and fixes
+# set boundary in main script because ffield is periodic
+units real
+# distribute electrode atoms among all processors:
+if "$(extract_setting(world_size) % 2) == 0" then "processors * * 2"
+
+atom_style full
+pair_style lj/cut/coul/long 12
+
+read_data "data.au-elyt"
+Reading data file ...
+  orthogonal box = (0 0 -10) to (1 1 10)
+  1 by 1 by 1 MPI processor grid
+  reading atoms ...
+  4 atoms
+Finding 1-2 1-3 1-4 neighbors ...
+  special bond factors lj:    0        0        0       
+  special bond factors coul:  0        0        0       
+     0 = max # of 1-2 neighbors
+     0 = max # of 1-3 neighbors
+     0 = max # of 1-4 neighbors
+     1 = max # of special neighbors
+  special bonds CPU = 0.000 seconds
+  read_data CPU = 0.003 seconds
+
+group bot type 1
+1 atoms in group bot
+group top type 2
+1 atoms in group top
+
+# get electrode charges
+variable q atom q
+compute qbot bot reduce sum v_q
+compute qtop top reduce sum v_q
+
+compute compute_pe all pe
+variable vpe equal c_compute_pe
+variable charge equal c_qtop
+fix fxprint all print 1 "${vpe}, ${charge}" file "out.csv"
+
+thermo_style custom step pe c_qbot c_qtop
+fix feta all property/atom d_eta ghost on
+set group bot d_eta 2.0
+Setting atom values ...
+  1 settings made for d_eta
+set group top d_eta 2.0
+Setting atom values ...
+  1 settings made for d_eta
+fix conp bot electrode/conp 0 2 couple top 1 symm on eta d_eta write_inv inv.csv write_vec vec.csv
+2 atoms in group conp_group
+
+run 0
+
+CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
+
+Your simulation uses code contributions which should be cited:
+
+- fix electrode command:
+
+@article{Ahrens2022
+author = {Ahrens-Iwers, Ludwig J.V. and Janssen, Mahijs and Tee, Shern R. and Mei{\ss}ner, Robert H.},
+doi = {10.1063/5.0099239},
+title = {{ELECTRODE: An electrochemistry package for LAMMPS}},
+journal = {The Journal of Chemical Physics},
+year = {2022}
+volume = {157},
+pages = {084801},
+}
+CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
+
+WARNING: No fixes with time integration, atoms won't move (src/verlet.cpp:60)
+Ewald/electrode initialization ...
+  using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
+WARNING: For better accuracy use 'pair_modify table 0' (src/kspace.cpp:365)
+  G vector (1/distance) = 0.32261103
+  estimated absolute RMS force accuracy = 3.8272011e-06
+  estimated relative force accuracy = 1.1525502e-08
+  KSpace vectors: actual max1d max3d = 52 50 515150
+                  kxmax kymax kzmax  = 1 1 50
+Generated 3 of 3 mixed pair_coeff terms from geometric mixing rule
+Neighbor list info ...
+  update: every = 1 steps, delay = 0 steps, check = yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 14
+  ghost atom cutoff = 14
+  binsize = 7, bins = 1 1 3
+  2 neighbor lists, perpetual/occasional/extra = 2 0 0
+  (1) pair lj/cut/coul/long, perpetual
+      attributes: half, newton on
+      pair build: half/bin/newton
+      stencil: half/bin/3d
+      bin: standard
+  (2) fix electrode/conp, perpetual, copy from (1)
+      attributes: half, newton on
+      pair build: copy
+      stencil: none
+      bin: none
+WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (src/domain.cpp:965)
+139.943964815502, 0.279214485147238
+Per MPI rank memory allocation (min/avg/max) = 144.2 | 144.2 | 144.2 Mbytes
+   Step         PotEng         c_qbot         c_qtop    
+         0   139.94396     -0.27921449     0.27921449   
+Loop time of 2.191e-06 on 1 procs for 0 steps with 4 atoms
+
+91.3% CPU use with 1 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 0          | 0          | 0          |   0.0 |  0.00
+Bond    | 0          | 0          | 0          |   0.0 |  0.00
+Kspace  | 0          | 0          | 0          |   0.0 |  0.00
+Neigh   | 0          | 0          | 0          |   0.0 |  0.00
+Comm    | 0          | 0          | 0          |   0.0 |  0.00
+Output  | 0          | 0          | 0          |   0.0 |  0.00
+Modify  | 0          | 0          | 0          |   0.0 |  0.00
+Other   |            | 2.191e-06  |            |       |100.00
+
+Nlocal:              4 ave           4 max           4 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:           3596 ave        3596 max        3596 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:           4790 ave        4790 max        4790 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 4790
+Ave neighs/atom = 1197.5
+Ave special neighs/atom = 0
+Neighbor list builds = 0
+Dangerous builds = 0
+
+Total wall time: 0:00:00
diff --git a/examples/PACKAGES/electrode/madelung/log.19Feb2024.eta_cg.g++.1 b/examples/PACKAGES/electrode/madelung/log.19Feb2024.eta_cg.g++.1
new file mode 100644
index 0000000000..edb2e434e6
--- /dev/null
+++ b/examples/PACKAGES/electrode/madelung/log.19Feb2024.eta_cg.g++.1
@@ -0,0 +1,139 @@
+LAMMPS (21 Nov 2023 - Development - patch_21Nov2023-668-g5b6c0c6b56)
+  using 1 OpenMP thread(s) per MPI task
+boundary p p f
+kspace_style ewald/electrode 1.0e-8
+kspace_modify slab 8.0 # ew3dc
+
+include "settings.mod" # styles, computes, groups and fixes
+# set boundary in main script because ffield is periodic
+units real
+# distribute electrode atoms among all processors:
+if "$(extract_setting(world_size) % 2) == 0" then "processors * * 2"
+
+atom_style full
+pair_style lj/cut/coul/long 12
+
+read_data "data.au-elyt"
+Reading data file ...
+  orthogonal box = (0 0 -10) to (1 1 10)
+  1 by 1 by 1 MPI processor grid
+  reading atoms ...
+  4 atoms
+Finding 1-2 1-3 1-4 neighbors ...
+  special bond factors lj:    0        0        0       
+  special bond factors coul:  0        0        0       
+     0 = max # of 1-2 neighbors
+     0 = max # of 1-3 neighbors
+     0 = max # of 1-4 neighbors
+     1 = max # of special neighbors
+  special bonds CPU = 0.000 seconds
+  read_data CPU = 0.003 seconds
+
+group bot type 1
+1 atoms in group bot
+group top type 2
+1 atoms in group top
+
+# get electrode charges
+variable q atom q
+compute qbot bot reduce sum v_q
+compute qtop top reduce sum v_q
+
+compute compute_pe all pe
+variable vpe equal c_compute_pe
+variable charge equal c_qtop
+fix fxprint all print 1 "${vpe}, ${charge}" file "out.csv"
+
+thermo_style custom step pe c_qbot c_qtop
+fix feta all property/atom d_eta ghost on
+set group bot d_eta 0.5
+Setting atom values ...
+  1 settings made for d_eta
+set group top d_eta 3.0
+Setting atom values ...
+  1 settings made for d_eta
+fix conp bot electrode/conp 0 2 couple top 1 symm on eta d_eta algo cg 1e-6
+2 atoms in group conp_group
+
+run 0
+
+CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
+
+Your simulation uses code contributions which should be cited:
+
+- fix electrode command:
+
+@article{Ahrens2022
+author = {Ahrens-Iwers, Ludwig J.V. and Janssen, Mahijs and Tee, Shern R. and Mei{\ss}ner, Robert H.},
+doi = {10.1063/5.0099239},
+title = {{ELECTRODE: An electrochemistry package for LAMMPS}},
+journal = {The Journal of Chemical Physics},
+year = {2022}
+volume = {157},
+pages = {084801},
+}
+CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
+
+WARNING: No fixes with time integration, atoms won't move (src/verlet.cpp:60)
+Ewald/electrode initialization ...
+  using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
+WARNING: For better accuracy use 'pair_modify table 0' (src/kspace.cpp:365)
+  G vector (1/distance) = 0.32261103
+  estimated absolute RMS force accuracy = 3.8272011e-06
+  estimated relative force accuracy = 1.1525502e-08
+  KSpace vectors: actual max1d max3d = 52 50 515150
+                  kxmax kymax kzmax  = 1 1 50
+Generated 3 of 3 mixed pair_coeff terms from geometric mixing rule
+Neighbor list info ...
+  update: every = 1 steps, delay = 0 steps, check = yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 14
+  ghost atom cutoff = 14
+  binsize = 7, bins = 1 1 3
+  2 neighbor lists, perpetual/occasional/extra = 2 0 0
+  (1) pair lj/cut/coul/long, perpetual
+      attributes: half, newton on
+      pair build: half/bin/newton
+      stencil: half/bin/3d
+      bin: standard
+  (2) fix electrode/conp, perpetual, copy from (1)
+      attributes: half, newton on
+      pair build: copy
+      stencil: none
+      bin: none
+WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (src/domain.cpp:965)
+165.519373910316, 0.29521534552818
+Per MPI rank memory allocation (min/avg/max) = 144.2 | 144.2 | 144.2 Mbytes
+   Step         PotEng         c_qbot         c_qtop    
+         0   165.51937     -0.29521535     0.29521535   
+Loop time of 2.797e-06 on 1 procs for 0 steps with 4 atoms
+
+71.5% CPU use with 1 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 0          | 0          | 0          |   0.0 |  0.00
+Bond    | 0          | 0          | 0          |   0.0 |  0.00
+Kspace  | 0          | 0          | 0          |   0.0 |  0.00
+Neigh   | 0          | 0          | 0          |   0.0 |  0.00
+Comm    | 0          | 0          | 0          |   0.0 |  0.00
+Output  | 0          | 0          | 0          |   0.0 |  0.00
+Modify  | 0          | 0          | 0          |   0.0 |  0.00
+Other   |            | 2.797e-06  |            |       |100.00
+
+Nlocal:              4 ave           4 max           4 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:           3596 ave        3596 max        3596 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:           4790 ave        4790 max        4790 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 4790
+Ave neighs/atom = 1197.5
+Ave special neighs/atom = 0
+Neighbor list builds = 0
+Dangerous builds = 0
+
+Average conjugate gradient steps: 1
+Total wall time: 0:00:00
diff --git a/examples/PACKAGES/electrode/madelung/log.19Feb2024.eta_mix.g++.1 b/examples/PACKAGES/electrode/madelung/log.19Feb2024.eta_mix.g++.1
new file mode 100644
index 0000000000..51eda0d870
--- /dev/null
+++ b/examples/PACKAGES/electrode/madelung/log.19Feb2024.eta_mix.g++.1
@@ -0,0 +1,138 @@
+LAMMPS (21 Nov 2023 - Development - patch_21Nov2023-668-g5b6c0c6b56)
+  using 1 OpenMP thread(s) per MPI task
+boundary p p f
+kspace_style ewald/electrode 1.0e-8
+kspace_modify slab 8.0 # ew3dc
+
+include "settings.mod" # styles, computes, groups and fixes
+# set boundary in main script because ffield is periodic
+units real
+# distribute electrode atoms among all processors:
+if "$(extract_setting(world_size) % 2) == 0" then "processors * * 2"
+
+atom_style full
+pair_style lj/cut/coul/long 12
+
+read_data "data.au-elyt"
+Reading data file ...
+  orthogonal box = (0 0 -10) to (1 1 10)
+  1 by 1 by 1 MPI processor grid
+  reading atoms ...
+  4 atoms
+Finding 1-2 1-3 1-4 neighbors ...
+  special bond factors lj:    0        0        0       
+  special bond factors coul:  0        0        0       
+     0 = max # of 1-2 neighbors
+     0 = max # of 1-3 neighbors
+     0 = max # of 1-4 neighbors
+     1 = max # of special neighbors
+  special bonds CPU = 0.000 seconds
+  read_data CPU = 0.003 seconds
+
+group bot type 1
+1 atoms in group bot
+group top type 2
+1 atoms in group top
+
+# get electrode charges
+variable q atom q
+compute qbot bot reduce sum v_q
+compute qtop top reduce sum v_q
+
+compute compute_pe all pe
+variable vpe equal c_compute_pe
+variable charge equal c_qtop
+fix fxprint all print 1 "${vpe}, ${charge}" file "out.csv"
+
+thermo_style custom step pe c_qbot c_qtop
+fix feta all property/atom d_eta ghost on
+set group bot d_eta 0.5
+Setting atom values ...
+  1 settings made for d_eta
+set group top d_eta 3.0
+Setting atom values ...
+  1 settings made for d_eta
+fix conp bot electrode/conp 0 2 couple top 1 symm on eta d_eta write_inv inv.csv write_vec vec.csv
+2 atoms in group conp_group
+
+run 0
+
+CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
+
+Your simulation uses code contributions which should be cited:
+
+- fix electrode command:
+
+@article{Ahrens2022
+author = {Ahrens-Iwers, Ludwig J.V. and Janssen, Mahijs and Tee, Shern R. and Mei{\ss}ner, Robert H.},
+doi = {10.1063/5.0099239},
+title = {{ELECTRODE: An electrochemistry package for LAMMPS}},
+journal = {The Journal of Chemical Physics},
+year = {2022}
+volume = {157},
+pages = {084801},
+}
+CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
+
+WARNING: No fixes with time integration, atoms won't move (src/verlet.cpp:60)
+Ewald/electrode initialization ...
+  using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
+WARNING: For better accuracy use 'pair_modify table 0' (src/kspace.cpp:365)
+  G vector (1/distance) = 0.32261103
+  estimated absolute RMS force accuracy = 3.8272011e-06
+  estimated relative force accuracy = 1.1525502e-08
+  KSpace vectors: actual max1d max3d = 52 50 515150
+                  kxmax kymax kzmax  = 1 1 50
+Generated 3 of 3 mixed pair_coeff terms from geometric mixing rule
+Neighbor list info ...
+  update: every = 1 steps, delay = 0 steps, check = yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 14
+  ghost atom cutoff = 14
+  binsize = 7, bins = 1 1 3
+  2 neighbor lists, perpetual/occasional/extra = 2 0 0
+  (1) pair lj/cut/coul/long, perpetual
+      attributes: half, newton on
+      pair build: half/bin/newton
+      stencil: half/bin/3d
+      bin: standard
+  (2) fix electrode/conp, perpetual, copy from (1)
+      attributes: half, newton on
+      pair build: copy
+      stencil: none
+      bin: none
+WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (src/domain.cpp:965)
+165.519373910316, 0.295215345528172
+Per MPI rank memory allocation (min/avg/max) = 144.2 | 144.2 | 144.2 Mbytes
+   Step         PotEng         c_qbot         c_qtop    
+         0   165.51937     -0.29521535     0.29521535   
+Loop time of 2.18e-06 on 1 procs for 0 steps with 4 atoms
+
+91.7% CPU use with 1 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 0          | 0          | 0          |   0.0 |  0.00
+Bond    | 0          | 0          | 0          |   0.0 |  0.00
+Kspace  | 0          | 0          | 0          |   0.0 |  0.00
+Neigh   | 0          | 0          | 0          |   0.0 |  0.00
+Comm    | 0          | 0          | 0          |   0.0 |  0.00
+Output  | 0          | 0          | 0          |   0.0 |  0.00
+Modify  | 0          | 0          | 0          |   0.0 |  0.00
+Other   |            | 2.18e-06   |            |       |100.00
+
+Nlocal:              4 ave           4 max           4 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:           3596 ave        3596 max        3596 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:           4790 ave        4790 max        4790 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 4790
+Ave neighs/atom = 1197.5
+Ave special neighs/atom = 0
+Neighbor list builds = 0
+Dangerous builds = 0
+
+Total wall time: 0:00:00
diff --git a/examples/PACKAGES/electrode/madelung/plate_cap.py b/examples/PACKAGES/electrode/madelung/plate_cap.py
index 62d52fe102..fcca166869 100755
--- a/examples/PACKAGES/electrode/madelung/plate_cap.py
+++ b/examples/PACKAGES/electrode/madelung/plate_cap.py
@@ -3,7 +3,6 @@
 import numpy as np
 from scipy.special import erf
 
-ETA = 2
 SQRT2 = np.sqrt(2)
 COULOMB = 332.06371  #  Coulomb constant in Lammps 'real' units
 QE2F = 23.060549
@@ -17,14 +16,14 @@ def lattice(length):
     return np.array(np.meshgrid(x, y)).T.reshape(-1, 2)
 
 
-def a_element(r):
+def a_element(r, eta):
     """Coulomb contribution of two Gaussians"""
-    return erf(ETA / SQRT2 * r) / r
+    return erf(eta * r) / r
 
 
-def b_element(r, q):
+def b_element(r, q, eta):
     """Coulomb contribution of a Gaussian with a point charge"""
-    return q * erf(ETA * r) / r
+    return q * erf(eta * r) / r
 
 
 a = 1  # nearest neighbor distance i.e. lattice constant / sqrt(2)
@@ -36,59 +35,65 @@ v = np.array([-0.5, 0.5]) * (QE2F / COULOMB)
 
 # distances to images within electrode and to opposite electrode
 distances = a * np.linalg.norm(lattice(LENGTH), axis=1)
-opposite_distances = np.sqrt(np.square(distances) + distance_plates ** 2)
+opposite_distances = np.sqrt(np.square(distances) + distance_plates**2)
 
-# self interaction and within original box
-A_11 = np.sqrt(2 / np.pi) * ETA
-A_12 = erf(ETA * distance_plates / SQRT2) / distance_plates
+for name, eta_elec in [("", [2.0, 2.0]), ("_eta_mix", [0.5, 3.0])]:
+    eta_mix = np.prod(eta_elec) / np.sqrt(np.sum(np.square(eta_elec)))
+    # self interaction and within original box
+    A_11 = np.sqrt(2 / np.pi) * eta_elec[0]
+    A_22 = np.sqrt(2 / np.pi) * eta_elec[1]
+    A_12 = erf(eta_mix * distance_plates) / distance_plates
 
-# interaction with periodic images
-A_11 += 4 * np.sum(a_element(distances))
-A_12 += 4 * np.sum(a_element(opposite_distances))
-A = np.array([[A_11, A_12], [A_12, A_11]])
-inv = np.linalg.inv(A)
-e = np.array([1, 1])
-inv -= np.matmul(inv, np.matmul(np.outer(e, e), inv)) / np.dot(e, np.dot(inv, e))
+    # interaction with periodic images
+    A_11 += 4 * np.sum(a_element(distances, eta_elec[0] / SQRT2))
+    A_22 += 4 * np.sum(a_element(distances, eta_elec[1] / SQRT2))
+    A_12 += 4 * np.sum(a_element(opposite_distances, eta_mix))
+    A = np.array([[A_11, A_12], [A_12, A_22]])
+    inv = np.linalg.inv(A)
+    e = np.array([1, 1])
+    inv -= np.matmul(inv, np.matmul(np.outer(e, e), inv)) / np.dot(e, np.dot(inv, e))
 
-# electrode-electrolyte interaction
-b = []
-for x in x_elec:
-    bi = 0
-    for y, q in zip(x_elyt, q_elyt):
-        d = abs(y - x)
-        bi += b_element(d, q)
-        image_distances = np.sqrt(np.square(distances) + d ** 2)
-        bi += 4 * np.sum(b_element(image_distances, q))
-    b.append(bi)
-b = np.array(b)
+    # electrode-electrolyte interaction
+    b = []
+    for x, eta in zip(x_elec, eta_elec):
+        bi = 0
+        for y, q in zip(x_elyt, q_elyt):
+            d = abs(y - x)
+            bi += b_element(d, q, eta)
+            image_distances = np.sqrt(np.square(distances) + d**2)
+            bi += 4 * np.sum(b_element(image_distances, q, eta))
+        b.append(bi)
+    b = np.array(b)
 
-# electrolyte-electrolyte energy
-elyt_11 = 4 * np.sum(1 / distances)
-distance_elyt = x_elyt[1] - x_elyt[0]
-elyt_12 = 1 / distance_elyt + 4 * np.sum(
-    1 / np.sqrt(np.square(distances) + distance_elyt ** 2)
-)
-elyt = np.array([[elyt_11, elyt_12], [elyt_12, elyt_11]])
-energy_elyt = 0.5 * np.dot(q_elyt, np.dot(elyt, q_elyt))
-
-# electrode charges and energy
-q = np.dot(inv, v - b)
-energy = COULOMB * (0.5 * np.dot(q, np.dot(A, q)) + np.dot(b, q) + energy_elyt)
-
-print(
-    "length, energy / kcal/mol, q1 / e, q2 / e, inv11 / A, inv12 / A, b1 / e/A, b2 / e/A"
-)
-print(
-    ", ".join(
-        [
-            str(LENGTH),
-            f"{energy:.8f}",
-            f"{q[0]:.10f}",
-            f"{q[1]:.10f}",
-            f"{inv[0, 0]:.10f}",
-            f"{inv[0, 1]:.10f}",
-            f"{b[0]:.8f}",
-            f"{b[1]:.8f}",
-        ]
+    # electrolyte-electrolyte energy
+    elyt_11 = 4 * np.sum(1 / distances)
+    distance_elyt = x_elyt[1] - x_elyt[0]
+    elyt_12 = 1 / distance_elyt + 4 * np.sum(
+        1 / np.sqrt(np.square(distances) + distance_elyt**2)
     )
-)
+    elyt = np.array([[elyt_11, elyt_12], [elyt_12, elyt_11]])
+    energy_elyt = 0.5 * np.dot(q_elyt, np.dot(elyt, q_elyt))
+
+    # electrode charges and energy
+    q = np.dot(inv, v - b)
+    energy = COULOMB * (0.5 * np.dot(q, np.dot(A, q)) + np.dot(b, q) + energy_elyt)
+
+    with open(f"plate_cap{name}.csv", "w") as f:
+        f.write(
+            "length, energy / kcal/mol, q1 / e, q2 / e, inv11 / A, inv12 / A, b1 / e/A, b2 / e/A\n"
+        )
+        f.write(
+            ", ".join(
+                [
+                    str(LENGTH),
+                    f"{energy:.8f}",
+                    f"{q[0]:.10f}",
+                    f"{q[1]:.10f}",
+                    f"{inv[0, 0]:.10f}",
+                    f"{inv[0, 1]:.10f}",
+                    f"{b[0]:.8f}",
+                    f"{b[1]:.8f}",
+                ]
+            )
+            + "\n"
+        )
diff --git a/examples/PACKAGES/electrode/madelung/test.sh b/examples/PACKAGES/electrode/madelung/test.sh
index edac04f5b1..a558ee6711 100644
--- a/examples/PACKAGES/electrode/madelung/test.sh
+++ b/examples/PACKAGES/electrode/madelung/test.sh
@@ -7,17 +7,27 @@ if [ ! -f $lmpbin ]; then
 fi
 
 ref_out="plate_cap.csv"
-if [ ! -f $ref_out ]; then
+ref_mix_out="plate_cap_eta_mix.csv"
+if [ ! -f $ref_out ] || [ ! -f $ref_mix_out ]; then
     echo "Generating reference data"
-    python3 plate_cap.py > $ref_out
+    python3 plate_cap.py
 fi
 
 echo "Running Lammps inputs"
+# w/o eta mixing
 rm -rf madelung.txt && touch madelung.txt
-for file in in.*; do
+for file in in.eta in.ewald-ew3dc in.ewald-ew2d in.pppm-ew3dc in.cg; do
     printf "\n$file\n" >> madelung.txt
     rm -f out.csv inv.csv vec.csv 
     $lmpbin -i $file &> /dev/null
     python3 eval.py $ref_out out.csv inv.csv vec.csv
 done
+
+# with eta mixing
+for file in in.eta_mix in.eta_cg; do
+    printf "\n$file\n" >> madelung.txt
+    rm -f out.csv inv.csv vec.csv 
+    $lmpbin -i $file &> /dev/null
+    python3 eval.py $ref_mix_out out.csv inv.csv vec.csv
+done
 cat madelung.txt
diff --git a/src/ELECTRODE/README b/src/ELECTRODE/README
new file mode 100644
index 0000000000..72a95b7fe9
--- /dev/null
+++ b/src/ELECTRODE/README
@@ -0,0 +1,17 @@
+This package provides the "fix electrode/*" commands which can be used in a
+LAMMPS input script. These fixes implement the constant potential method, which
+minimizes the energy of electrodes as a function of atom charges at given
+electric potentials or electrode charges.
+
+See the doc page for the fix electrode/conp command to get started.  There are
+example scripts for using this package in examples/PACKAGES/electrode.
+
+This package uses an external library in lib/electrode which must be compiled
+before making LAMMPS. For a CMake build, the location of the LAPACK library
+should be linked automatically. Alternatively, the "USE_INTERNAL_LINALG" option
+can be used to enable the bundled library. See the doc page on "Packages with
+extra build options" for more information.
+
+The primary people who created this package are Ludwig Ahrens-Iwers, Shern Tee
+(s.tee@griffith.edu.au) and Robert Meißner (robert.meissner@tuhh.de). Contact
+them directly if you have questions.
diff --git a/src/ELECTRODE/electrode_math.h b/src/ELECTRODE/electrode_math.h
index 5992df2289..4a3cb7bac4 100644
--- a/src/ELECTRODE/electrode_math.h
+++ b/src/ELECTRODE/electrode_math.h
@@ -18,22 +18,20 @@
 #ifndef LMP_ELECTRODE_MATH_H
 #define LMP_ELECTRODE_MATH_H
 
+#include "ewald_const.h"
 #include "math_const.h"
 
+#include <cmath>
+
 namespace LAMMPS_NS {
+using namespace EwaldConst;
 
 namespace ElectrodeMath {
-  static constexpr double EWALD_P = 0.3275911;
-  static constexpr double A1 = 0.254829592;
-  static constexpr double A2 = -0.284496736;
-  static constexpr double A3 = 1.421413741;
-  static constexpr double A4 = -1.453152027;
-  static constexpr double A5 = 1.061405429;
   static constexpr double ERFCMAX = 5.8;    // erfc(ERFCMAX) < machine epsilon(double)
 
   static double safe_erfc(double x)
   {
-    if (x > ERFCMAX) return 0.;
+    if (x > ERFCMAX) return 0.0;
     double expm2 = exp(-x * x);
     double t = 1.0 / (1.0 + EWALD_P * x);
     return t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2;
@@ -42,14 +40,14 @@ namespace ElectrodeMath {
   static double safe_derfcr(double x, double &erfc)
   {
     if (x > ERFCMAX) {
-      erfc = 0.;
-      return 0.;
+      erfc = 0.0;
+      return 0.0;
     }
     double x2 = x * x;
     double expm2 = exp(-x2);
     double t = 1.0 / (1.0 + EWALD_P * x);
     erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2;
-    return -erfc - 2 * expm2 * x / MathConst::MY_PIS;
+    return -erfc - 2.0 * expm2 * x / MathConst::MY_PIS;
   }
 }    // namespace ElectrodeMath
 
diff --git a/src/ELECTRODE/electrode_matrix.cpp b/src/ELECTRODE/electrode_matrix.cpp
index 86761742d4..9c2da1d13b 100644
--- a/src/ELECTRODE/electrode_matrix.cpp
+++ b/src/ELECTRODE/electrode_matrix.cpp
@@ -43,6 +43,7 @@ ElectrodeMatrix::ElectrodeMatrix(LAMMPS *lmp, int electrode_group, double eta) :
   groupbit = group->bitmask[igroup];
   ngroup = group->count(igroup);
   this->eta = eta;
+  etaflag = false;
   tfflag = false;
 }
 
@@ -72,6 +73,14 @@ void ElectrodeMatrix::setup_tf(const std::map<int, double> &tf_types)
 
 /* ---------------------------------------------------------------------- */
 
+void ElectrodeMatrix::setup_eta(int index)
+{
+  etaflag = true;
+  eta_index = index;
+}
+
+/* ---------------------------------------------------------------------- */
+
 void ElectrodeMatrix::compute_array(double **array, bool timer_flag)
 {
   // setting all entries of coulomb matrix to zero
@@ -115,8 +124,6 @@ void ElectrodeMatrix::pair_contribution(double **array)
   int nlocal = atom->nlocal;
   int newton_pair = force->newton_pair;
 
-  double const etaij = eta * eta / sqrt(2.0 * eta * eta);    // see mw ewald theory eq. (29)-(30)
-
   // neighbor list will be ready because called from post_neighbor
   inum = list->inum;
   ilist = list->ilist;
@@ -135,6 +142,7 @@ void ElectrodeMatrix::pair_contribution(double **array)
     xtmp = x[i][0];
     ytmp = x[i][1];
     ztmp = x[i][2];
+    double const eta_i = etaflag ? atom->dvector[eta_index][i] : eta;
     itype = type[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
@@ -152,6 +160,9 @@ void ElectrodeMatrix::pair_contribution(double **array)
       jtype = type[j];
 
       if (rsq < cutsq[itype][jtype]) {
+        double const eta_j = etaflag ? atom->dvector[eta_index][j] : eta;
+        double const etaij = eta_i * eta_j / sqrt(eta_i * eta_i + eta_j * eta_j);
+
         r = sqrt(rsq);
         rinv = 1.0 / r;
         aij = rinv;
@@ -178,7 +189,10 @@ void ElectrodeMatrix::self_contribution(double **array)
   const double preta = MY_SQRT2 / MY_PIS;
 
   for (int i = 0; i < nlocal; i++)
-    if (mask[i] & groupbit) { array[mpos[i]][mpos[i]] += preta * eta - selfint; }
+    if (mask[i] & groupbit) {
+      double const eta_i = etaflag ? atom->dvector[eta_index][i] : eta;
+      array[mpos[i]][mpos[i]] += preta * eta_i - selfint;
+    }
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/ELECTRODE/electrode_matrix.h b/src/ELECTRODE/electrode_matrix.h
index 8499cfdb34..1c64d8a4c4 100644
--- a/src/ELECTRODE/electrode_matrix.h
+++ b/src/ELECTRODE/electrode_matrix.h
@@ -30,6 +30,7 @@ class ElectrodeMatrix : protected Pointers {
   ElectrodeMatrix(class LAMMPS *, int, double);
   void setup(const std::unordered_map<tagint, int> &, class Pair *, class NeighList *);
   void setup_tf(const std::map<int, double> &);
+  void setup_eta(int);
   void compute_array(double **, bool);
   int igroup;
 
@@ -39,6 +40,8 @@ class ElectrodeMatrix : protected Pointers {
   double **cutsq;
   double g_ewald, eta;
   bool tfflag;
+  bool etaflag;
+  int eta_index;
   std::map<int, double> tf_types;
   std::unordered_map<tagint, int> tag_to_iele;
   bool assigned;
diff --git a/src/ELECTRODE/electrode_vector.cpp b/src/ELECTRODE/electrode_vector.cpp
index 8511ddc17c..fc2cca5e46 100644
--- a/src/ELECTRODE/electrode_vector.cpp
+++ b/src/ELECTRODE/electrode_vector.cpp
@@ -29,6 +29,7 @@
 #include "neigh_list.h"
 #include "pair.h"
 
+#include <cassert>
 #include <cmath>
 #include <exception>
 
@@ -47,6 +48,7 @@ ElectrodeVector::ElectrodeVector(LAMMPS *lmp, int sensor_group, int source_group
   source_grpbit = group->bitmask[source_group];
   this->eta = eta;
   tfflag = false;
+  etaflag = false;
 
   kspace_time_total = 0;
   pair_time_total = 0;
@@ -93,6 +95,14 @@ void ElectrodeVector::setup_tf(const std::map<int, double> &tf_types)
 
 /* ---------------------------------------------------------------------- */
 
+void ElectrodeVector::setup_eta(int index)
+{
+  etaflag = true;
+  eta_index = index;
+}
+
+/* ---------------------------------------------------------------------- */
+
 void ElectrodeVector::compute_vector(double *vector)
 {
   MPI_Barrier(world);
@@ -121,7 +131,6 @@ void ElectrodeVector::compute_vector(double *vector)
 
 void ElectrodeVector::pair_contribution(double *vector)
 {
-  double const etaij = eta * MY_ISQRT2;
   double **x = atom->x;
   double *q = atom->q;
   int *type = atom->type;
@@ -142,6 +151,7 @@ void ElectrodeVector::pair_contribution(double *vector)
     double const xtmp = x[i][0];
     double const ytmp = x[i][1];
     double const ztmp = x[i][2];
+    double const eta_i = etaflag ? atom->dvector[eta_index][i] : eta;
     int itype = type[i];
     int *jlist = firstneigh[i];
     int jnum = numneigh[i];
@@ -158,18 +168,22 @@ void ElectrodeVector::pair_contribution(double *vector)
       double const rsq = delx * delx + dely * dely + delz * delz;
       int jtype = type[j];
       if (rsq >= cutsq[itype][jtype]) continue;
+      double const eta_j = etaflag ? atom->dvector[eta_index][j] : eta;
+      double etaij;
+      if (i_in_sensor && j_in_sensor) {
+        etaij = eta_i * eta_j / sqrt(eta_i * eta_i + eta_j * eta_j);
+      } else if (i_in_sensor) {
+        etaij = eta_i;
+      } else {
+        assert(j_in_sensor);
+        etaij = eta_j;
+      }
       double const r = sqrt(rsq);
       double const rinv = 1.0 / r;
       double aij = rinv;
       aij *= ElectrodeMath::safe_erfc(g_ewald * r);
-      if (invert_source)
-        aij -= ElectrodeMath::safe_erfc(eta * r) * rinv;
-      else
-        aij -= ElectrodeMath::safe_erfc(etaij * r) * rinv;
-      if (i_in_sensor) {
-        vector[i] += aij * q[j];
-        //} else if (j_in_sensor) {
-      }
+      aij -= ElectrodeMath::safe_erfc(etaij * r) * rinv;
+      if (i_in_sensor) { vector[i] += aij * q[j]; }
       if (j_in_sensor && (!invert_source || !i_in_sensor)) { vector[j] += aij * q[i]; }
     }
   }
@@ -189,9 +203,10 @@ void ElectrodeVector::self_contribution(double *vector)
 
   for (int ii = 0; ii < inum; ii++) {
     int const i = ilist[ii];
+    double const eta_i = etaflag ? atom->dvector[eta_index][i] : eta;
     bool const i_in_sensor = (mask[i] & groupbit);
     bool const i_in_source = !!(mask[i] & source_grpbit) != invert_source;
-    if (i_in_sensor && i_in_source) vector[i] += (preta * eta - selfint) * q[i];
+    if (i_in_sensor && i_in_source) vector[i] += (preta * eta_i - selfint) * q[i];
   }
 }
 
diff --git a/src/ELECTRODE/electrode_vector.h b/src/ELECTRODE/electrode_vector.h
index e7f637dd2d..a4f274a049 100644
--- a/src/ELECTRODE/electrode_vector.h
+++ b/src/ELECTRODE/electrode_vector.h
@@ -29,6 +29,7 @@ class ElectrodeVector : protected Pointers {
   ~ElectrodeVector() override;
   void setup(class Pair *, class NeighList *, bool);
   void setup_tf(const std::map<int, double> &);
+  void setup_eta(int);
   void compute_vector(double *);
   int igroup, source_group;
 
@@ -39,6 +40,8 @@ class ElectrodeVector : protected Pointers {
   double **cutsq;
   double g_ewald, eta;
   bool tfflag;
+  bool etaflag;
+  int eta_index;
   std::map<int, double> tf_types;
   class Pair *pair;
   class NeighList *list;
diff --git a/src/ELECTRODE/fix_electrode_conp.cpp b/src/ELECTRODE/fix_electrode_conp.cpp
index e8b11c330b..94c085de5c 100644
--- a/src/ELECTRODE/fix_electrode_conp.cpp
+++ b/src/ELECTRODE/fix_electrode_conp.cpp
@@ -98,26 +98,30 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) :
   top_group = 0;
   intelflag = false;
   tfflag = false;
+  etaflag = false;
   timer_flag = false;
 
   update_time = 0;
   mult_time = 0;
   n_call = n_cg_step = 0;
 
+  qtotal = 0.;
+  qtotal_var_style = VarStyle::UNSET;
+
   // read fix command
   fixname = std::string(arg[0]);
   groups = std::vector<int>(1, igroup);
   group_bits = std::vector<int>(1, groupbit);
   group_psi_var_names = std::vector<std::string>(1);
   group_psi_var_styles = std::vector<VarStyle>(1, VarStyle::CONST);
-  group_psi = std::vector<double>(1);
+  group_psi_const = std::vector<double>(1);
   etypes_neighlists = false;
   if (strstr(arg[3], "v_") == arg[3]) {
     std::string vname = arg[3];
     group_psi_var_names[0] = vname.substr(2);
     group_psi_var_styles[0] = VarStyle::EQUAL;
   } else
-    group_psi[0] = utils::numeric(FLERR, arg[3], false, lmp);
+    group_psi_const[0] = utils::numeric(FLERR, arg[3], false, lmp);
   char *eta_str = arg[4];
   eta = utils::numeric(FLERR, eta_str, false, lmp);
   int iarg = 5;
@@ -133,12 +137,12 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) :
         std::string vname = arg[iarg];
         group_psi_var_names.push_back(vname.substr(2));
         group_psi_var_styles.push_back(VarStyle::EQUAL);
-        group_psi.push_back(0.);
+        group_psi_const.push_back(0.);
       } else {
         std::string null;
         group_psi_var_names.push_back(null);
         group_psi_var_styles.push_back(VarStyle::CONST);
-        group_psi.push_back(utils::numeric(FLERR, arg[iarg], false, lmp));
+        group_psi_const.push_back(utils::numeric(FLERR, arg[iarg], false, lmp));
       }
     } else if ((strcmp(arg[iarg], "algo") == 0)) {
       if (!default_algo) error->one(FLERR, "Algorithm can be set only once");
@@ -196,8 +200,32 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) :
       thermo_temp = force->boltz / force->qe2f * utils::numeric(FLERR, arg[++iarg], false, lmp);
       thermo_time = utils::numeric(FLERR, arg[++iarg], false, lmp);
       thermo_init = utils::inumeric(FLERR, arg[++iarg], false, lmp);
-      // toggle parameters
-    } else if ((strcmp(arg[iarg], "etypes") == 0)) {
+    } else if ((strcmp(arg[iarg], "qtotal") == 0)) {
+      if (iarg + 2 > narg) error->all(FLERR, "Need one argument after qtotal keyword");
+      if (strcmp(this->style, "electrode/conq") == 0)
+        error->all(FLERR, "qtotal keyword not available for electrode/conq");
+      ++iarg;
+      if (strstr(arg[iarg], "v_") == arg[iarg]) {
+        std::string vname = arg[iarg];
+        qtotal_var_name = vname.substr(2);
+        qtotal_var_style = VarStyle::EQUAL;
+      } else {
+        qtotal = utils::numeric(FLERR, arg[iarg], false, lmp);
+        qtotal_var_style = VarStyle::CONST;
+      }
+    } else if ((strcmp(arg[iarg], "eta") == 0)) {
+      if (iarg + 2 > narg) error->all(FLERR, "Need two arguments after eta command");
+      etaflag = true;
+      int is_double, cols, ghost;
+      eta_index = atom->find_custom_ghost(arg[++iarg] + 2, is_double, cols, ghost);
+      if (eta_index == -1)
+        error->all(FLERR, "eta keyword requires name of previously defined property");
+      if (!is_double) error->all(FLERR, "eta keyword requires double-valued property/atom vector");
+      if (cols != 0) error->all(FLERR, "eta keyword requires property/atom vector not an array");
+      if (!ghost) error->all(FLERR, "eta keyword requires property/atom fix with ghost on");
+    }
+    // toggle parameters
+    else if ((strcmp(arg[iarg], "etypes") == 0)) {
       etypes_neighlists = utils::logical(FLERR, arg[++iarg], false, lmp);
     } else if ((strncmp(arg[iarg], "symm", 4) == 0)) {
       symm = utils::logical(FLERR, arg[++iarg], false, lmp);
@@ -209,6 +237,12 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) :
     iarg++;
   }
 
+  if (qtotal_var_style != VarStyle::UNSET) {
+    if (symm) error->all(FLERR, "{} cannot use qtotal keyword with symm on", this->style);
+  }
+
+  // computatonal potential
+  group_psi = std::vector<double>(groups.size());
   // union of all coupled groups
   std::string union_group = "conp_group";
   std::string group_cmd = union_group + " union";
@@ -226,6 +260,7 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) :
   if (need_elec_vector) elec_vector = new ElectrodeVector(lmp, igroup, igroup, eta, false);
   assert(groups.size() == group_bits.size());
   assert(groups.size() == group_psi.size());
+  assert(groups.size() == group_psi_const.size());
   assert(groups.size() == group_psi_var_styles.size());
   assert(groups.size() == group_psi_var_names.size());
   assert(igroup == elyt_vector->igroup);
@@ -375,6 +410,31 @@ void FixElectrodeConp::init()
     if (strncmp(modify->fix[i]->style, "electrode", 9) == 0) count++;
   if (count > 1) error->all(FLERR, "More than one fix electrode");
 
+  // make sure electrode atoms are not integrated if a matrix is used for electrode-electrode interaction
+  int const nlocal = atom->nlocal;
+  int *mask = atom->mask;
+  Fix **fix = modify->fix;
+  if (matrix_algo) {
+    std::vector<char *> integrate_ids = std::vector<char *>();
+    for (int i = 0; i < modify->nfix; i++) {
+      if (fix[i]->time_integrate == 0) continue;
+      int electrode_mover = 0;
+      int fix_groupbit = fix[i]->groupbit;
+      for (int j = 0; j < nlocal; j++)
+        if ((mask[j] & fix_groupbit) && (mask[j] & groupbit)) electrode_mover = 1;
+      MPI_Allreduce(MPI_IN_PLACE, &electrode_mover, 1, MPI_INT, MPI_SUM, world);
+      if (electrode_mover && comm->me == 0) integrate_ids.push_back(fix[i]->id);
+    }
+    if (comm->me == 0)
+      for (char *fix_id : integrate_ids)
+        error->warning(
+            FLERR,
+            "Electrode atoms are integrated by fix {}, but fix electrode is using a matrix method. "
+            "For "
+            "mobile electrodes use the conjugate gradient algorithm without matrix ('algo cg').",
+            fix_id);
+  }
+
   // check for package intel
   if (etypes_neighlists)
     request_etypes_neighlists();
@@ -449,6 +509,7 @@ void FixElectrodeConp::setup_post_neighbor()
   // get equal-style variable ids:
   group_psi_var_ids = std::vector<int>(num_of_groups, -1);
   for (int g = 0; g < num_of_groups; g++) {
+    assert(group_psi_var_styles[g] != VarStyle::UNSET);
     if (group_psi_var_styles[g] == VarStyle::CONST) continue;
     const char *var_name = group_psi_var_names[g].c_str();
     int var_id = input->variable->find(var_name);
@@ -457,13 +518,23 @@ void FixElectrodeConp::setup_post_neighbor()
       error->all(FLERR, "Variable '{}' for fix {} is not equal-style", var_name, style);
     group_psi_var_ids[g] = var_id;
   }
+  if (qtotal_var_style == VarStyle::EQUAL) {
+    const char *var_name = qtotal_var_name.c_str();
+    int var_id = input->variable->find(var_name);
+    if (var_id < 0) error->all(FLERR, "Variable '{}' for fix electrode does not exist", var_name);
+    if (!input->variable->equalstyle(var_id))
+      error->all(FLERR, "Variable '{}' for fix electrode is not equal-style", var_name);
+    qtotal_var_id = var_id;
+  }
 
   // pair and list setups:
 
   evscale = force->qe2f / force->qqrd2e;
   elyt_vector->setup(pair, vec_neighlist, timer_flag);
+  if (etaflag) elyt_vector->setup_eta(eta_index);
   if (need_elec_vector) {
     elec_vector->setup(pair, mat_neighlist, timer_flag);
+    if (etaflag) elec_vector->setup_eta(eta_index);
     if (tfflag) elec_vector->setup_tf(tf_types);
   }
 
@@ -498,7 +569,8 @@ void FixElectrodeConp::setup_post_neighbor()
       if (etypes_neighlists) neighbor->build_one(mat_neighlist, 0);
       auto array_compute = std::unique_ptr<ElectrodeMatrix>(new ElectrodeMatrix(lmp, igroup, eta));
       array_compute->setup(tag_to_iele, pair, mat_neighlist);
-      if (tfflag) { array_compute->setup_tf(tf_types); }
+      if (etaflag) array_compute->setup_eta(eta_index);
+      if (tfflag) array_compute->setup_tf(tf_types);
       array_compute->compute_array(elastance, timer_flag);
     }    // write_mat before proceeding
     if (comm->me == 0 && write_mat) {
@@ -805,6 +877,8 @@ void FixElectrodeConp::update_charges()
     }
     MPI_Allreduce(MPI_IN_PLACE, &sb_charges.front(), num_of_groups, MPI_DOUBLE, MPI_SUM, world);
     update_psi();    // use for equal-style and conq
+    if (qtotal_var_style != VarStyle::UNSET)
+      update_psi_qtotal();    // use for qtotal; same for thermo
     for (int g = 0; g < num_of_groups; g++)
       for (int j = 0; j < nlocalele; j++) q_local[j] += sd_vectors[g][list_iele[j]] * group_psi[g];
     MPI_Barrier(world);
@@ -907,12 +981,22 @@ std::vector<double> FixElectrodeConp::gather_ngroup(std::vector<double> x_local)
 }
 
 /* ----------------------------------------------------------------------
-   ensure total electrode charge is 0 if symm
+   ensure total electrode charge is 0 if symm and qtotal if qtotal is used
 ------------------------------------------------------------------------- */
 
 std::vector<double> FixElectrodeConp::constraint_correction(std::vector<double> x)
 {
-  return constraint_projection(std::move(x));
+  if (symm || qtotal_var_style != VarStyle::UNSET) {
+    if (qtotal_var_style == VarStyle::EQUAL) qtotal = input->variable->compute_equal(qtotal_var_id);
+    double sum = 0.;
+    for (double xi : x) sum += xi;
+    MPI_Allreduce(MPI_IN_PLACE, &sum, 1, MPI_DOUBLE, MPI_SUM, world);
+    if (qtotal_var_style != VarStyle::UNSET) sum -= qtotal;
+    sum /= ngroup;
+    for (double &xi : x) xi -= sum;
+    return x;
+  }
+  return x;
 }
 
 /* ----------------------------------------------------------------------
@@ -921,7 +1005,7 @@ std::vector<double> FixElectrodeConp::constraint_correction(std::vector<double>
 
 std::vector<double> FixElectrodeConp::constraint_projection(std::vector<double> x)
 {
-  if (symm) {
+  if (symm || qtotal_var_style != VarStyle::UNSET) {
     double sum = 0.;
     for (double xi : x) sum += xi;
     MPI_Allreduce(MPI_IN_PLACE, &sum, 1, MPI_DOUBLE, MPI_SUM, world);
@@ -979,13 +1063,28 @@ std::vector<double> FixElectrodeConp::times_elastance(std::vector<double> x)
 void FixElectrodeConp::update_psi()
 {
   for (int g = 0; g < num_of_groups; g++) {
-    if (group_psi_var_styles[g] == VarStyle::CONST) continue;
-    group_psi[g] = input->variable->compute_equal(group_psi_var_ids[g]);
+    if (group_psi_var_styles[g] == VarStyle::CONST)
+      group_psi[g] = group_psi_const[g];
+    else
+      group_psi[g] = input->variable->compute_equal(group_psi_var_ids[g]);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
+void FixElectrodeConp::update_psi_qtotal()
+{
+  if (qtotal_var_style == VarStyle::EQUAL) qtotal = input->variable->compute_equal(qtotal_var_id);
+  double q_current = 0.;
+  for (int i = 0; i < num_of_groups; i++) {
+    q_current += sb_charges[i];
+    for (int j = 0; j < num_of_groups; j++) q_current += macro_capacitance[i][j] * group_psi[j];
+  }
+  double add_psi = (qtotal - q_current) / macro_capacitance_sum;
+  for (int i = 0; i < num_of_groups; i++) group_psi[i] += add_psi;
+}
+/* ---------------------------------------------------------------------- */
+
 void FixElectrodeConp::compute_macro_matrices()
 {
   assert(algo == Algo::MATRIX_INV);
@@ -1041,6 +1140,10 @@ void FixElectrodeConp::compute_macro_matrices()
       }
     }
   }
+
+  macro_capacitance_sum = 0.;
+  for (int i = 0; i < num_of_groups; i++)
+    for (int j = 0; j < num_of_groups; j++) macro_capacitance_sum += macro_capacitance[i][j];
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1097,7 +1200,7 @@ double FixElectrodeConp::self_energy(int eflag)
   // corrections to energy due to self interaction
   double const qqrd2e = force->qqrd2e;
   int const nlocal = atom->nlocal;
-  double const pre = eta / sqrt(MY_2PI) * qqrd2e;
+  double const pre = 1. / sqrt(MY_2PI) * qqrd2e;
   int *mask = atom->mask;
   int *type = atom->type;
   double *q = atom->q;
@@ -1105,7 +1208,8 @@ double FixElectrodeConp::self_energy(int eflag)
   for (int i = 0; i < nlocal; i++) {
     if (groupbit & mask[i]) {
       double const q2 = q[i] * q[i];
-      double e = pre * q2;
+      double ieta = etaflag ? atom->dvector[eta_index][i] : eta;
+      double e = ieta * pre * q2;
       if (tfflag && (groupbit & mask[i])) e += 0.5 * qqrd2e * q2 * tf_types[type[i]];
       energy += e;
       if (eflag) {
@@ -1145,6 +1249,7 @@ double FixElectrodeConp::gausscorr(int eflag, bool fflag)
     double xtmp = x[i][0];
     double ytmp = x[i][1];
     double ztmp = x[i][2];
+    double const eta_i = etaflag ? atom->dvector[eta_index][i] : eta;
     int itype = type[i];
     int *jlist = firstneigh[i];
     int jnum = numneigh[i];
@@ -1153,7 +1258,6 @@ double FixElectrodeConp::gausscorr(int eflag, bool fflag)
       int const j = jlist[jj] & NEIGHMASK;
       bool j_in_ele = groupbit & mask[j];
       if (!(i_in_ele || j_in_ele)) continue;
-      double eta_ij = (i_in_ele && j_in_ele) ? eta / MY_SQRT2 : eta;
 
       double delx = xtmp - x[j][0];
       double dely = ytmp - x[j][1];
@@ -1162,6 +1266,16 @@ double FixElectrodeConp::gausscorr(int eflag, bool fflag)
       int jtype = type[j];
 
       if (rsq < force->pair->cutsq[itype][jtype]) {
+        double const eta_j = etaflag ? atom->dvector[eta_index][j] : eta;
+        double eta_ij;
+        if (i_in_ele && j_in_ele)
+          eta_ij = eta_i * eta_j / sqrt(eta_i * eta_i + eta_j * eta_j);
+        else if (i_in_ele)
+          eta_ij = eta_i;
+        else {
+          assert(j_in_ele);
+          eta_ij = eta_j;
+        }
         double r2inv = 1.0 / rsq;
         double r = sqrt(rsq);
         double erfc_etar = 0.;
diff --git a/src/ELECTRODE/fix_electrode_conp.h b/src/ELECTRODE/fix_electrode_conp.h
index 46da1eec35..a1d7530bd1 100644
--- a/src/ELECTRODE/fix_electrode_conp.h
+++ b/src/ELECTRODE/fix_electrode_conp.h
@@ -70,7 +70,7 @@ class FixElectrodeConp : public Fix {
 
  protected:
   enum class Algo { MATRIX_INV, MATRIX_CG, CG };
-  enum class VarStyle { CONST, EQUAL };
+  enum class VarStyle { CONST, EQUAL, UNSET };
   virtual void update_psi();
   virtual void pre_update(){};
   virtual void recompute_potential(std::vector<double>, std::vector<double>){};
@@ -79,6 +79,7 @@ class FixElectrodeConp : public Fix {
   virtual void compute_macro_matrices();
   std::vector<double> ele_ele_interaction(const std::vector<double> &);
   std::vector<double> group_psi;
+  std::vector<double> group_psi_const;    // needed to undo qtotal psi updates
   std::vector<int> group_bits;
   std::vector<int> groups;
   int num_of_groups;
@@ -100,6 +101,10 @@ class FixElectrodeConp : public Fix {
   std::string fixname;    // used by electrode/ffield to set up internal efield
   bool intelflag;
   inline virtual void intel_pack_buffers() {}
+  double qtotal;
+  std::string qtotal_var_name;
+  int qtotal_var_id;
+  VarStyle qtotal_var_style;
 
  private:
   std::string output_file_inv, output_file_mat, output_file_vec;
@@ -132,6 +137,8 @@ class FixElectrodeConp : public Fix {
   int get_top_group();    // used by ffield
   int top_group;          // used by ffield
   bool tfflag;
+  int eta_index;    // index of atom property for eta
+  bool etaflag;     // eta specified as atom property
   bool timer_flag;
   std::map<int, double> tf_types;
   // cg
@@ -142,6 +149,9 @@ class FixElectrodeConp : public Fix {
   std::vector<double> gather_ngroup(std::vector<double>);
   std::vector<double> gather_elevec_local(ElectrodeVector *);
   void set_charges(std::vector<double>);
+  // qtotal
+  double macro_capacitance_sum;
+  void update_psi_qtotal();
 
   // fix-specific electrode ID storage system:
 
diff --git a/src/ELECTRODE/fix_electrode_conq.cpp b/src/ELECTRODE/fix_electrode_conq.cpp
index 0d3d1d2aaf..a6baa1e122 100644
--- a/src/ELECTRODE/fix_electrode_conq.cpp
+++ b/src/ELECTRODE/fix_electrode_conq.cpp
@@ -30,7 +30,7 @@ FixElectrodeConq::FixElectrodeConq(LAMMPS *lmp, int narg, char **arg) :
     FixElectrodeConp(lmp, narg, arg)
 {
   // copy const-style values across because update_psi will change group_psi
-  group_q = group_psi;
+  group_q = group_psi_const;
 
   if (symm) {
     if (num_of_groups == 1)
diff --git a/src/ELECTRODE/fix_electrode_thermo.cpp b/src/ELECTRODE/fix_electrode_thermo.cpp
index 343bf14069..92db4b3ee0 100644
--- a/src/ELECTRODE/fix_electrode_thermo.cpp
+++ b/src/ELECTRODE/fix_electrode_thermo.cpp
@@ -47,7 +47,8 @@ FixElectrodeThermo::FixElectrodeThermo(LAMMPS *lmp, int narg, char **arg) :
   if (thermo_time < SMALL) error->all(FLERR, "Keyword temp not set or zero in electrode/thermo");
 
   thermo_random = new RanMars(lmp, thermo_init);
-  if (group_psi_var_styles[0] == VarStyle::CONST) delta_psi_0 = group_psi[1] - group_psi[0];
+  if (group_psi_var_styles[0] == VarStyle::CONST)
+    delta_psi_0 = group_psi_const[1] - group_psi_const[0];
 }
 
 /* ----------------------------------------------------------------------- */
@@ -102,7 +103,7 @@ void FixElectrodeThermo::update_psi()
   double const delta_psi = group_psi_old[1] - group_psi_old[0];
 
   // target potential difference from input parameters
-  if (group_psi_var_styles[0] != VarStyle::CONST) {
+  if (group_psi_var_styles[0] == VarStyle::EQUAL) {
     delta_psi_0 = input->variable->compute_equal(group_psi_var_ids[1]) -
         input->variable->compute_equal(group_psi_var_ids[0]);
   }
diff --git a/src/ELECTRODE/pppm_electrode.cpp b/src/ELECTRODE/pppm_electrode.cpp
index 98a76a7cca..ee34def74d 100644
--- a/src/ELECTRODE/pppm_electrode.cpp
+++ b/src/ELECTRODE/pppm_electrode.cpp
@@ -438,6 +438,7 @@ void PPPMElectrode::compute(int eflag, int vflag)
 
   start_compute();
 
+  /*
   if (compute_vector_called && last_invert_source) {
     // electrolyte_density_brick is filled, so we can grab only electrode atoms.
     // Does not work for direct cg algorithm because electrode charges change after compute_vector.
@@ -453,15 +454,17 @@ void PPPMElectrode::compute(int eflag, int vflag)
           density_brick[nz][ny][nx] += electrolyte_density_brick[nz][ny][nx];
         }
   } else {
-    make_rho();
+  */
+  particle_map();
+  make_rho();
 
-    // all procs communicate density values from their ghost cells
-    //   to fully sum contribution in their 3d bricks
-    // remap from 3d decomposition to FFT decomposition
+  // all procs communicate density values from their ghost cells
+  //   to fully sum contribution in their 3d bricks
+  // remap from 3d decomposition to FFT decomposition
 
-    gc->reverse_comm(Grid3d::KSPACE, this, REVERSE_RHO, 1, sizeof(FFT_SCALAR), gc_buf1, gc_buf2,
-                     MPI_FFT_SCALAR);
-  }
+  gc->reverse_comm(Grid3d::KSPACE, this, REVERSE_RHO, 1, sizeof(FFT_SCALAR), gc_buf1, gc_buf2,
+                   MPI_FFT_SCALAR);
+  //}
 
   brick2fft();
 
@@ -584,6 +587,7 @@ void PPPMElectrode::compute_vector(double *vec, int sensor_grpbit, int source_gr
   // electrolyte density (without writing an additional function)
   FFT_SCALAR ***density_brick_real = density_brick;
   FFT_SCALAR *density_fft_real = density_fft;
+  particle_map();
   make_rho_in_brick(source_grpbit, electrolyte_density_brick, invert_source);
   density_brick = electrolyte_density_brick;
   density_fft = electrolyte_density_fft;
@@ -669,7 +673,8 @@ void PPPMElectrode::compute_matrix(bigint *imat, double **matrix, bool timer_fla
   // fft green's function k -> r (double)
   double *greens_real;
   memory->create(greens_real, nz_pppm * ny_pppm * nx_pppm, "pppm/electrode:greens_real");
-  memset(greens_real, 0, (std::size_t)nz_pppm * (std::size_t)ny_pppm * (std::size_t)nx_pppm * sizeof(double));
+  memset(greens_real, 0,
+         (std::size_t) nz_pppm * (std::size_t) ny_pppm * (std::size_t) nx_pppm * sizeof(double));
   for (int i = 0, n = 0; i < nfft; i++) {
     work2[n++] = greensfn[i];
     work2[n++] = ZEROF;
@@ -862,7 +867,7 @@ void PPPMElectrode::two_step_multiplication(bigint *imat, double *greens_real, d
 
   double **gw;
   memory->create(gw, nmat, nxyz, "pppm/electrode:gw");
-  memset(&(gw[0][0]), 0, (std::size_t)nmat * (std::size_t)nxyz * sizeof(double));
+  memset(&(gw[0][0]), 0, (std::size_t) nmat * (std::size_t) nxyz * sizeof(double));
 
   auto fmod = [](int x, int n) {    // fast unsigned mod
     int r = abs(x);
@@ -981,17 +986,18 @@ void PPPMElectrode::allocate()
   // returns local owned and ghost grid bounds
   // setup communication patterns and buffers
 
-  gc = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm,
-                  nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
-                  nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
+  gc = new Grid3d(lmp, world, nx_pppm, ny_pppm, nz_pppm, nxlo_in, nxhi_in, nylo_in, nyhi_in,
+                  nzlo_in, nzhi_in, nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out, nzhi_out);
 
-  gc->setup_comm(ngc_buf1,ngc_buf2);
+  gc->setup_comm(ngc_buf1, ngc_buf2);
 
-  if (differentiation_flag) npergrid = 1;
-  else npergrid = 3;
+  if (differentiation_flag)
+    npergrid = 1;
+  else
+    npergrid = 3;
 
-  memory->create(gc_buf1,npergrid*ngc_buf1,"pppm:gc_buf1");
-  memory->create(gc_buf2,npergrid*ngc_buf2,"pppm:gc_buf2");
+  memory->create(gc_buf1, npergrid * ngc_buf1, "pppm:gc_buf1");
+  memory->create(gc_buf2, npergrid * ngc_buf2, "pppm:gc_buf2");
 
   // tally local grid sizes
   // ngrid = count of owned+ghost grid cells on this proc
@@ -1000,67 +1006,63 @@ void PPPMElectrode::allocate()
   // nfft = FFT points in x-pencil FFT decomposition on this proc
   // nfft_both = greater of nfft and nfft_brick
 
-  ngrid = (nxhi_out-nxlo_out+1) * (nyhi_out-nylo_out+1) *
-    (nzhi_out-nzlo_out+1);
+  ngrid = (nxhi_out - nxlo_out + 1) * (nyhi_out - nylo_out + 1) * (nzhi_out - nzlo_out + 1);
 
-  nfft_brick = (nxhi_in-nxlo_in+1) * (nyhi_in-nylo_in+1) *
-    (nzhi_in-nzlo_in+1);
+  nfft_brick = (nxhi_in - nxlo_in + 1) * (nyhi_in - nylo_in + 1) * (nzhi_in - nzlo_in + 1);
 
-  nfft = (nxhi_fft-nxlo_fft+1) * (nyhi_fft-nylo_fft+1) *
-    (nzhi_fft-nzlo_fft+1);
+  nfft = (nxhi_fft - nxlo_fft + 1) * (nyhi_fft - nylo_fft + 1) * (nzhi_fft - nzlo_fft + 1);
 
-  nfft_both = MAX(nfft,nfft_brick);
+  nfft_both = MAX(nfft, nfft_brick);
 
   // allocate distributed grid data
 
-  memory->create3d_offset(density_brick,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                          nxlo_out,nxhi_out,"pppm:density_brick");
+  memory->create3d_offset(density_brick, nzlo_out, nzhi_out, nylo_out, nyhi_out, nxlo_out, nxhi_out,
+                          "pppm:density_brick");
 
-  memory->create(density_fft,nfft_both,"pppm:density_fft");
-  memory->create(greensfn,nfft_both,"pppm:greensfn");
-  memory->create(work1,2*nfft_both,"pppm:work1");
-  memory->create(work2,2*nfft_both,"pppm:work2");
-  memory->create(vg,nfft_both,6,"pppm:vg");
+  memory->create(density_fft, nfft_both, "pppm:density_fft");
+  memory->create(greensfn, nfft_both, "pppm:greensfn");
+  memory->create(work1, 2 * nfft_both, "pppm:work1");
+  memory->create(work2, 2 * nfft_both, "pppm:work2");
+  memory->create(vg, nfft_both, 6, "pppm:vg");
 
   if (triclinic == 0) {
-    memory->create1d_offset(fkx,nxlo_fft,nxhi_fft,"pppm:fkx");
-    memory->create1d_offset(fky,nylo_fft,nyhi_fft,"pppm:fky");
-    memory->create1d_offset(fkz,nzlo_fft,nzhi_fft,"pppm:fkz");
+    memory->create1d_offset(fkx, nxlo_fft, nxhi_fft, "pppm:fkx");
+    memory->create1d_offset(fky, nylo_fft, nyhi_fft, "pppm:fky");
+    memory->create1d_offset(fkz, nzlo_fft, nzhi_fft, "pppm:fkz");
   } else {
-    memory->create(fkx,nfft_both,"pppm:fkx");
-    memory->create(fky,nfft_both,"pppm:fky");
-    memory->create(fkz,nfft_both,"pppm:fkz");
+    memory->create(fkx, nfft_both, "pppm:fkx");
+    memory->create(fky, nfft_both, "pppm:fky");
+    memory->create(fkz, nfft_both, "pppm:fkz");
   }
 
   if (differentiation_flag == 1) {
-    memory->create3d_offset(u_brick,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                          nxlo_out,nxhi_out,"pppm:u_brick");
+    memory->create3d_offset(u_brick, nzlo_out, nzhi_out, nylo_out, nyhi_out, nxlo_out, nxhi_out,
+                            "pppm:u_brick");
 
-    memory->create(sf_precoeff1,nfft_both,"pppm:sf_precoeff1");
-    memory->create(sf_precoeff2,nfft_both,"pppm:sf_precoeff2");
-    memory->create(sf_precoeff3,nfft_both,"pppm:sf_precoeff3");
-    memory->create(sf_precoeff4,nfft_both,"pppm:sf_precoeff4");
-    memory->create(sf_precoeff5,nfft_both,"pppm:sf_precoeff5");
-    memory->create(sf_precoeff6,nfft_both,"pppm:sf_precoeff6");
+    memory->create(sf_precoeff1, nfft_both, "pppm:sf_precoeff1");
+    memory->create(sf_precoeff2, nfft_both, "pppm:sf_precoeff2");
+    memory->create(sf_precoeff3, nfft_both, "pppm:sf_precoeff3");
+    memory->create(sf_precoeff4, nfft_both, "pppm:sf_precoeff4");
+    memory->create(sf_precoeff5, nfft_both, "pppm:sf_precoeff5");
+    memory->create(sf_precoeff6, nfft_both, "pppm:sf_precoeff6");
 
   } else {
-    memory->create3d_offset(vdx_brick,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                            nxlo_out,nxhi_out,"pppm:vdx_brick");
-    memory->create3d_offset(vdy_brick,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                            nxlo_out,nxhi_out,"pppm:vdy_brick");
-    memory->create3d_offset(vdz_brick,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                            nxlo_out,nxhi_out,"pppm:vdz_brick");
+    memory->create3d_offset(vdx_brick, nzlo_out, nzhi_out, nylo_out, nyhi_out, nxlo_out, nxhi_out,
+                            "pppm:vdx_brick");
+    memory->create3d_offset(vdy_brick, nzlo_out, nzhi_out, nylo_out, nyhi_out, nxlo_out, nxhi_out,
+                            "pppm:vdy_brick");
+    memory->create3d_offset(vdz_brick, nzlo_out, nzhi_out, nylo_out, nyhi_out, nxlo_out, nxhi_out,
+                            "pppm:vdz_brick");
   }
 
   // summation coeffs
 
   order_allocated = order;
-  if (!stagger_flag) memory->create(gf_b,order,"pppm:gf_b");
-  memory->create2d_offset(rho1d,3,-order/2,order/2,"pppm:rho1d");
-  memory->create2d_offset(drho1d,3,-order/2,order/2,"pppm:drho1d");
-  memory->create2d_offset(rho_coeff,order,(1-order)/2,order/2,"pppm:rho_coeff");
-  memory->create2d_offset(drho_coeff,order,(1-order)/2,order/2,
-                          "pppm:drho_coeff");
+  if (!stagger_flag) memory->create(gf_b, order, "pppm:gf_b");
+  memory->create2d_offset(rho1d, 3, -order / 2, order / 2, "pppm:rho1d");
+  memory->create2d_offset(drho1d, 3, -order / 2, order / 2, "pppm:drho1d");
+  memory->create2d_offset(rho_coeff, order, (1 - order) / 2, order / 2, "pppm:rho_coeff");
+  memory->create2d_offset(drho_coeff, order, (1 - order) / 2, order / 2, "pppm:drho_coeff");
 
   // create 2 FFTs and a Remap
   // 1st FFT keeps data in FFT decomposition
@@ -1069,20 +1071,17 @@ void PPPMElectrode::allocate()
 
   int tmp;
 
-  fft1 = new FFT3d(lmp,world,nx_pppm,ny_pppm,nz_pppm,
-                   nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
-                   nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
-                   0,0,&tmp,collective_flag);
+  fft1 = new FFT3d(lmp, world, nx_pppm, ny_pppm, nz_pppm, nxlo_fft, nxhi_fft, nylo_fft, nyhi_fft,
+                   nzlo_fft, nzhi_fft, nxlo_fft, nxhi_fft, nylo_fft, nyhi_fft, nzlo_fft, nzhi_fft,
+                   0, 0, &tmp, collective_flag);
 
-  fft2 = new FFT3d(lmp,world,nx_pppm,ny_pppm,nz_pppm,
-                   nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
-                   nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
-                   0,0,&tmp,collective_flag);
+  fft2 = new FFT3d(lmp, world, nx_pppm, ny_pppm, nz_pppm, nxlo_fft, nxhi_fft, nylo_fft, nyhi_fft,
+                   nzlo_fft, nzhi_fft, nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in, 0, 0,
+                   &tmp, collective_flag);
 
-  remap = new Remap(lmp,world,
-                    nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
-                    nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
-                    1,0,0,FFT_PRECISION,collective_flag);
+  remap = new Remap(lmp, world, nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in, nxlo_fft,
+                    nxhi_fft, nylo_fft, nyhi_fft, nzlo_fft, nzhi_fft, 1, 0, 0, FFT_PRECISION,
+                    collective_flag);
 
   // ELECTRODE specific allocations
 
diff --git a/unittest/force-styles/tests/in.conp b/unittest/force-styles/tests/in.conp
index 92d2f63cd1..08673ec20b 100644
--- a/unittest/force-styles/tests/in.conp
+++ b/unittest/force-styles/tests/in.conp
@@ -22,3 +22,4 @@ angle_coeff *
 
 group bot type 1
 group top type 2
+group ele type 1 2
diff --git a/unittest/force-styles/tests/kspace-ewald_conp_charge.yaml b/unittest/force-styles/tests/kspace-ewald_conp_charge.yaml
index d0eaee6b92..9bc190a766 100644
--- a/unittest/force-styles/tests/kspace-ewald_conp_charge.yaml
+++ b/unittest/force-styles/tests/kspace-ewald_conp_charge.yaml
@@ -1,6 +1,6 @@
 ---
-lammps_version: 23 Jun 2022
-date_generated: Wed Sep 21 13:52:53 2022
+lammps_version: 7 Feb 2024
+date_generated: Mon Mar  4 09:44:30 2024
 epsilon: 1e-12
 skip_tests: gpu kokkos_omp omp
 prerequisites: ! |
@@ -16,6 +16,7 @@ post_commands: ! |
   kspace_modify gewald 0.23118
   kspace_modify slab ew2d
   fix fxcpm bot electrode/conp -1.0 1.805 couple top 1.0 symm on
+  fix fxforce ele setforce 0 0 0
 input_file: in.conp
 pair_style: coul/long 15.0
 pair_coeff: ! |
@@ -27,97 +28,97 @@ init_coul: 2.215589572896434
 init_stress: ! |2-
    0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 init_forces: ! |2
-    1  2.0780648532795694e-04  1.9949672015209204e-03  3.1005914149473996e+00
-    2 -1.6777235182686288e-02  2.1481432256290419e-03  3.0881659196467988e+00
-    3  6.0082164895554737e-04  5.1573260226633801e-03  3.1029192412328555e+00
-    4 -1.6728974802490675e-02  6.1174723156886242e-03  3.0909324782862346e+00
-    5  4.2029366155132378e-02 -2.3455526736195693e-03 -1.5659617577954634e+00
-    6  5.5635790919204904e-02 -2.4542947062522369e-03 -1.5693827709331334e+00
-    7  4.2014920784252008e-02 -7.5287470219125008e-04 -1.5671265392163820e+00
-    8  5.5808767852333470e-02 -9.9105389808573120e-04 -1.5707104957299389e+00
-    9 -5.0959878750421551e-02 -2.3630298689785601e-03 -1.5769250181497101e+00
-   10 -3.3526564930579039e-02 -2.3802275431282884e-03 -1.5617801011657175e+00
-   11 -5.1236396351794389e-02 -4.9531100598979201e-04 -1.5779995894034005e+00
-   12 -3.3740693032952060e-02 -1.0210406243572182e-03 -1.5630986537874150e+00
-   13 -1.1437102611353016e-03 -4.6454866413029015e-05  5.4282837980149448e-03
-   14  2.3914999373115431e-03 -1.6478680244651469e-04  2.9802178734319239e-02
-   15  3.9287193302652786e-05 -2.5715673267285659e-05  2.8944525105129479e-03
-   16  2.0458480716482328e-03 -1.2119161321908735e-04  3.3689550843809452e-02
-   17 -2.7146073277767471e-03 -8.2376243258224663e-04  2.6564130941474612e-02
-   18  1.3669692885198135e-03 -4.2357196145489820e-04  3.2396141113926739e-02
-   19  3.0143371860819995e-04 -8.6218593339583785e-04  2.6284521141350669e-02
-   20  1.1542435168435056e-03 -2.7252318260838826e-04  3.4237916528138110e-02
-   21 -1.2350056952573553e-03  4.8655691135364269e-04  5.9284283442393631e-03
-   22  2.3656743884722890e-03  9.6575340844312705e-04  2.9811074931784823e-02
-   23  4.6754986244969657e-05  3.0149464050350903e-04  3.4630785686112129e-03
-   24  2.0301227080749633e-03  6.3879578068684812e-04  3.3653437189053413e-02
-   25 -2.3656211013513076e-03 -8.0454594828768334e-04  2.8476980555362911e-02
-   26  1.1566723797447039e-03 -3.9614599888570504e-04  3.2873323713155905e-02
-   27  2.8784994028036400e-04 -8.3661697184444898e-04  2.8317655886021253e-02
-   28  9.3882364605486020e-04 -2.3327601777843495e-04  3.4334676606415648e-02
-   29 -4.7969977052124917e-04 -1.2933334305373028e-04 -1.2336987392568071e-02
-   30  6.4733118786851766e-05 -1.3190918849005797e-04 -1.2737933567178844e-02
-   31  2.4269094157913586e-04 -1.3093943526788584e-04 -1.2136133260085013e-02
-   32  1.7452552740941527e-04 -1.1792779046242341e-04 -1.4181538324619835e-02
-   33 -3.8366266481516803e-04 -7.1061854758754556e-05 -1.3699106365426135e-02
-   34  2.8849004082563746e-05 -5.7838605310673531e-05 -1.3764181266896890e-02
-   35  2.2648059665862587e-04 -7.2851385190891320e-05 -1.3537361892926607e-02
-   36  1.2929221129083645e-04 -4.2862960950045859e-05 -1.4926105930886896e-02
-   37 -4.7698025941707008e-04  2.9971529466656788e-04 -1.2393604822896313e-02
-   38  6.4231095731188766e-05  2.7548977518460050e-04 -1.2789498345723021e-02
-   39  2.4169204779864826e-04  3.0552093685810269e-04 -1.2193908285665961e-02
-   40  1.7324998349441456e-04  2.2898000918153004e-04 -1.4225267020837207e-02
-   41 -3.4345772150395188e-04 -9.8012060153887415e-05 -1.4482722052972283e-02
-   42  2.0345466940577010e-05 -8.5250083485342566e-05 -1.4497101004472062e-02
-   43  2.0917627239292995e-04 -1.0051271468149899e-04 -1.4335313646556430e-02
-   44  1.1456796622437295e-04 -6.7553675788598551e-05 -1.5543196158604005e-02
+    1  2.0780648532797705e-04  1.9949672015210172e-03  3.1005914149473988e+00
+    2 -1.6777235182686239e-02  2.1481432256291962e-03  3.0881659196467979e+00
+    3  6.0082164895566240e-04  5.1573260226632518e-03  3.1029192412328550e+00
+    4 -1.6728974802490665e-02  6.1174723156885227e-03  3.0909324782862333e+00
+    5  4.2029366155132364e-02 -2.3455526736195940e-03 -1.5659617577954636e+00
+    6  5.5635790919204904e-02 -2.4542947062522933e-03 -1.5693827709331336e+00
+    7  4.2014920784251980e-02 -7.5287470219118655e-04 -1.5671265392163822e+00
+    8  5.5808767852333491e-02 -9.9105389808568046e-04 -1.5707104957299391e+00
+    9 -5.0959878750421538e-02 -2.3630298689785965e-03 -1.5769250181497101e+00
+   10 -3.3526564930579081e-02 -2.3802275431283552e-03 -1.5617801011657175e+00
+   11 -5.1236396351794437e-02 -4.9531100598972577e-04 -1.5779995894034005e+00
+   12 -3.3740693032952088e-02 -1.0210406243571744e-03 -1.5630986537874152e+00
+   13  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   14  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   15  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   16  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   17  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   18  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   19  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   20  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   21  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   22  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   23  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   24  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   25  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   26  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   27  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   28  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   29  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   30  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   31  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   32  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   33  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   34  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   35  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   36  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   37  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   38  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   39  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   40  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   41  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   42  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   43  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   44  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 run_vdwl: 0
-run_coul: 6.662694629990089
+run_coul: 6.662694556930397
 run_stress: ! |2-
    0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 run_forces: ! |2
-    1  2.4590612609445102e-04  1.9614041218568861e-03  3.0874291949281147e+00
-    2 -1.6646393233505193e-02  2.1136941574790400e-03  3.0751132981100078e+00
-    3  6.3535217476586373e-04  5.1012487117746350e-03  3.0897301317927290e+00
-    4 -1.6598677148313409e-02  6.0535154567108685e-03  3.0778495346381409e+00
-    5  4.2257888534896988e-02 -2.3028533365965051e-03 -1.5593100596807521e+00
-    6  5.5690508027606708e-02 -2.4104721624763235e-03 -1.5626897542440843e+00
-    7  4.2243973130370149e-02 -7.6153220413259775e-04 -1.5604618260035832e+00
-    8  5.5862703939049158e-02 -9.9825803703216718e-04 -1.5640031014325448e+00
-    9 -5.1059409954744304e-02 -2.3195553026588347e-03 -1.5701677851024036e+00
-   10 -3.3824298857146967e-02 -2.3375522139358631e-03 -1.5551647619109401e+00
-   11 -5.1334079184640377e-02 -5.0583705005136689e-04 -1.5712298444761112e+00
-   12 -3.4037363466305925e-02 -1.0275978089057873e-03 -1.5564691026885336e+00
-   13 -1.1767076011504501e-03 -4.7681963272732406e-05  5.6177800544716262e-03
-   14  2.3826294437743331e-03 -1.6349140124449633e-04  2.9836275824428962e-02
-   15  4.1635367214843796e-05 -2.7287283914685102e-05  3.0967012748694773e-03
-   16  2.0334138778396313e-03 -1.1986039772787527e-04  3.3679727731055195e-02
-   17 -2.7026084826568797e-03 -8.1815079485725360e-04  2.6574868248546435e-02
-   18  1.3603406762441243e-03 -4.1902595052936860e-04  3.2373613783594497e-02
-   19  2.9940467686436986e-04 -8.5646794759970863e-04  2.6296909514905095e-02
-   20  1.1469475577225402e-03 -2.6907093945665336e-04  3.4197589258157073e-02
-   21 -1.2662578686531134e-03  4.9770031968890231e-04  6.1141650873547037e-03
-   22  2.3568409810395020e-03  9.5671841592085381e-04  2.9844352492872490e-02
-   23  4.8965507374742117e-05  3.1616579858329929e-04  3.6613180489820005e-03
-   24  2.0177821554069170e-03  6.3083810187911310e-04  3.3643353017422439e-02
-   25 -2.3537455003017457e-03 -7.9846295760147956e-04  2.8468250829639500e-02
-   26  1.1507655048236000e-03 -3.9159985067612060e-04  3.2839870487003708e-02
-   27  2.8582562554448814e-04 -8.3038492818152999e-04  2.8309777443009273e-02
-   28  9.3274285761092680e-04 -2.2997823984283208e-04  3.4287630335266286e-02
-   29 -4.7502414048888327e-04 -1.2847214455389489e-04 -1.2453998829891042e-02
-   30  6.3675154563755000e-05 -1.3104204562344653e-04 -1.2848240218511071e-02
-   31  2.4068203429808906e-04 -1.3007692195448562e-04 -1.2254443488117142e-02
-   32  1.7286880375665112e-04 -1.1713944254614034e-04 -1.4279748536149278e-02
-   33 -3.7975304094097439e-04 -7.0481503989778179e-05 -1.3805675914786045e-02
-   34  2.8172580422574018e-05 -5.7414812953860394e-05 -1.3867744309035916e-02
-   35  2.2448042926710853e-04 -7.2258368796568220e-05 -1.3645037267085249e-02
-   36  1.2804302797390986e-04 -4.2547023166907131e-05 -1.5019384166440140e-02
-   37 -4.7231656130619122e-04  2.9759414220405656e-04 -1.2510308090752414e-02
-   38  6.3181133759042352e-05  2.7366531591578288e-04 -1.2899551129187765e-02
-   39  2.3968383978790597e-04  3.0338116638639894e-04 -1.2311919409216509e-02
-   40  1.7159804905743131e-04  2.2744508340904917e-04 -1.4323246783684936e-02
-   41 -3.3986046011853923e-04 -9.7274166717457145e-05 -1.4583591616566860e-02
-   42  1.9775995840422978e-05 -8.4675117269364898e-05 -1.4595729099929700e-02
-   43  2.0725822207231137e-04 -9.9764702521265135e-05 -1.4437236974797743e-02
-   44  1.1345006523014360e-04 -6.7103771021504086e-05 -1.5632251527466189e-02
+    1  2.4590683616588566e-04  1.9614031577864796e-03  3.0874292323214036e+00
+    2 -1.6646391056109287e-02  2.1136931488459688e-03  3.0751133344119794e+00
+    3  6.3535284261976926e-04  5.1012486499620825e-03  3.0897301693141319e+00
+    4 -1.6598674977089559e-02  6.0535154273402107e-03  3.0778495712965124e+00
+    5  4.2257890248609438e-02 -2.3028528771736243e-03 -1.5593100743454236e+00
+    6  5.5690508840794850e-02 -2.4104716736066006e-03 -1.5626897684026868e+00
+    7  4.2243974848895645e-02 -7.6153214117129057e-04 -1.5604618407191602e+00
+    8  5.5862704771247280e-02 -9.9825804058299890e-04 -1.5640031157582626e+00
+    9 -5.1059412079591346e-02 -2.3195548367148894e-03 -1.5701677996288432e+00
+   10 -3.3824301875619085e-02 -2.3375517362837230e-03 -1.5551647759943836e+00
+   11 -5.1334081278836916e-02 -5.0583707575582921e-04 -1.5712298590599834e+00
+   12 -3.4037366496628307e-02 -1.0275977984217173e-03 -1.5564691169384559e+00
+   13  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   14  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   15  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   16  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   17  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   18  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   19  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   20  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   21  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   22  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   23  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   24  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   25  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   26  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   27  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   28  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   29  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   30  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   31  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   32  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   33  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   34  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   35  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   36  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   37  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   38  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   39  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   40  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   41  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   42  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   43  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   44  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 ...
diff --git a/unittest/force-styles/tests/kspace-pppm_conp_charge.yaml b/unittest/force-styles/tests/kspace-pppm_conp_charge.yaml
index 0555d12ed4..38161ee044 100644
--- a/unittest/force-styles/tests/kspace-pppm_conp_charge.yaml
+++ b/unittest/force-styles/tests/kspace-pppm_conp_charge.yaml
@@ -1,6 +1,6 @@
 ---
-lammps_version: 23 Jun 2022
-date_generated: Wed Sep 21 13:52:39 2022
+lammps_version: 7 Feb 2024
+date_generated: Mon Mar  4 09:44:31 2024
 epsilon: 3e-12
 skip_tests: gpu kokkos_omp omp
 prerequisites: ! |
@@ -16,6 +16,7 @@ post_commands: ! |
   kspace_modify gewald 0.23118
   kspace_modify slab 3.0
   fix fxcpm bot electrode/conp -1.0 1.805 couple top 1.0 symm on
+  fix fxforce ele setforce 0 0 0
 input_file: in.conp
 pair_style: coul/long 15.0
 pair_coeff: ! |
@@ -23,101 +24,101 @@ pair_coeff: ! |
 extract: ! ""
 natoms: 44
 init_vdwl: 0
-init_coul: 2.2156402256727614
+init_coul: 2.215640225672775
 init_stress: ! |2-
    0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 init_forces: ! |2
-    1  2.0996096688279944e-04  1.9837586784580306e-03  3.1004822661058822e+00
-    2 -1.6783332510617883e-02  2.1368843599407611e-03  3.0880130470329230e+00
-    3  6.0300296042517466e-04  5.1688381279905342e-03  3.1028182137891114e+00
-    4 -1.6735061532950901e-02  6.1290626039690339e-03  3.0907879891042755e+00
-    5  4.2014131860757913e-02 -2.3478381081742388e-03 -1.5658874682481487e+00
-    6  5.5659823770659422e-02 -2.4566144388410410e-03 -1.5693278833316506e+00
-    7  4.1999624791768865e-02 -7.5066218795240259e-04 -1.5670569809441617e+00
-    8  5.5832732887661884e-02 -9.8883264742169940e-04 -1.5706605160409139e+00
-    9 -5.0976953599115804e-02 -2.3653810185280950e-03 -1.5768945194236066e+00
-   10 -3.3513771125456657e-02 -2.3824712764543426e-03 -1.5616806812004886e+00
-   11 -5.1253442064492741e-02 -4.9304425051529275e-04 -1.5779738349804424e+00
-   12 -3.3727836471637192e-02 -1.0188844490582761e-03 -1.5630041309277038e+00
-   13 -1.1453068449918257e-03 -4.7335833322794788e-05  5.4292779404649470e-03
-   14  2.3900993287279790e-03 -1.6878550058260119e-04  2.9808528147740175e-02
-   15  4.0078428215627730e-05 -2.6184607051201481e-05  2.8941780231019881e-03
-   16  2.0473353699190459e-03 -1.2552209515766962e-04  3.3684989110502959e-02
-   17 -2.7210216747431955e-03 -8.2349543008294359e-04  2.6567504438257068e-02
-   18  1.3656002828979516e-03 -4.2323438710338486e-04  3.2404938547366383e-02
-   19  3.0785575286292939e-04 -8.6186674263511191e-04  2.6288541663855129e-02
-   20  1.1555469330548321e-03 -2.7230960410720359e-04  3.4235148032534163e-02
-   21 -1.2368093613861506e-03  4.8760847861882366e-04  5.9296798954256557e-03
-   22  2.3643140421916085e-03  9.6975102599399746e-04  2.9817231402721564e-02
-   23  4.7705653522709085e-05  3.0203836842154655e-04  3.4631818106649337e-03
-   24  2.0316297431160258e-03  6.4335031755788927e-04  3.3648629802522749e-02
-   25 -2.3728144718995455e-03 -8.0497592536520339e-04  2.8474707915345274e-02
-   26  1.1555985481661916e-03 -3.9649433660109970e-04  3.2876098209196375e-02
-   27  2.9459292459149998e-04 -8.3700881746301306e-04  2.8316136079038545e-02
-   28  9.4027352090446912e-04 -2.3371025598546553e-04  3.4325153603153732e-02
-   29 -5.2133856931286127e-04 -1.4498587872629142e-04 -1.2345168780426297e-02
-   30  7.0344538924238829e-05 -1.4805754895657979e-04 -1.2765142487049358e-02
-   31  2.7857644686035687e-04 -1.4667349483298643e-04 -1.2140095836769501e-02
-   32  1.7479865631996218e-04 -1.3335074368636031e-04 -1.4152171307753206e-02
-   33 -4.2607366742485959e-04 -7.2661709209033136e-05 -1.3713642029394900e-02
-   34  3.4224570995904750e-05 -5.9352423088202727e-05 -1.3797063100154012e-02
-   35  2.6332638434852483e-04 -7.4573909050170201e-05 -1.3548481212572565e-02
-   36  1.2956589882656260e-04 -4.4469640559387641e-05 -1.4903750420442119e-02
-   37 -5.1855202137555164e-04  3.1768672664149645e-04 -1.2402413876588990e-02
-   38  6.9744637010106164e-05  2.9379923413403201e-04 -1.2817238930047764e-02
-   39  2.7753815671547851e-04  3.2364879063853462e-04 -1.2198687630220518e-02
-   40  1.7359211286247436e-04  2.4661766514410969e-04 -1.4196607980261094e-02
-   41 -3.8501793941197521e-04 -9.8818538537990245e-05 -1.4472863167957050e-02
-   42  2.5779944975443997e-05 -8.5975255466014692e-05 -1.4504933662725301e-02
-   43  2.4535595442148292e-04 -1.0127263490049206e-04 -1.4321758087878972e-02
-   44  1.1457678622968817e-04 -6.8200688092186871e-05 -1.5499407046729242e-02
+    1  2.0996096688288336e-04  1.9837586784579777e-03  3.1004822661058866e+00
+    2 -1.6783332510618088e-02  2.1368843599407052e-03  3.0880130470329270e+00
+    3  6.0300296042527755e-04  5.1688381279906869e-03  3.1028182137891140e+00
+    4 -1.6735061532951116e-02  6.1290626039692056e-03  3.0907879891042778e+00
+    5  4.2014131860757857e-02 -2.3478381081742123e-03 -1.5658874682481501e+00
+    6  5.5659823770659547e-02 -2.4566144388410201e-03 -1.5693278833316524e+00
+    7  4.1999624791768810e-02 -7.5066218795247393e-04 -1.5670569809441637e+00
+    8  5.5832732887661961e-02 -9.8883264742177161e-04 -1.5706605160409157e+00
+    9 -5.0976953599115846e-02 -2.3653810185280759e-03 -1.5768945194236084e+00
+   10 -3.3513771125456573e-02 -2.3824712764543218e-03 -1.5616806812004898e+00
+   11 -5.1253442064492769e-02 -4.9304425051536007e-04 -1.5779738349804435e+00
+   12 -3.3727836471637081e-02 -1.0188844490583557e-03 -1.5630041309277052e+00
+   13  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   14  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   15  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   16  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   17  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   18  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   19  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   20  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   21  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   22  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   23  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   24  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   25  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   26  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   27  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   28  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   29  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   30  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   31  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   32  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   33  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   34  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   35  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   36  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   37  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   38  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   39  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   40  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   41  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   42  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   43  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   44  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 run_vdwl: 0
-run_coul: 6.662844717848837
+run_coul: 6.662844644802024
 run_stress: ! |2-
    0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 run_forces: ! |2
-    1  2.4838374656870440e-04  1.9503798034564181e-03  3.0873204052231675e+00
-    2 -1.6652792550963507e-02  2.1026197438206527e-03  3.0749612313228378e+00
-    3  6.3785681925848106e-04  5.1125747842690368e-03  3.0896293583134606e+00
-    4 -1.6605065971975488e-02  6.0649203428150876e-03  3.0777057441315305e+00
-    5  4.2242720963296274e-02 -2.3050719334346786e-03 -1.5592361839228894e+00
-    6  5.5714272244614366e-02 -2.4127241272198356e-03 -1.5626351053810128e+00
-    7  4.2228744507856318e-02 -7.5938680980830215e-04 -1.5603926173148708e+00
-    8  5.5886400463161408e-02 -9.9610508394880446e-04 -1.5639532943562677e+00
-    9 -5.1076260774663269e-02 -2.3218376339398856e-03 -1.5701372580807029e+00
-   10 -3.3811558405906099e-02 -2.3397301017438945e-03 -1.5550660295253020e+00
-   11 -5.1350900750795785e-02 -5.0363945411724615e-04 -1.5712039970018119e+00
-   12 -3.4024561045577786e-02 -1.0255079139096881e-03 -1.5563752017113657e+00
-   13 -1.1783355827329363e-03 -4.8584425213699084e-05  5.6187810005819084e-03
-   14  2.3812433011153062e-03 -1.6745821939573657e-04  2.9842604329757397e-02
-   15  4.2472634867315832e-05 -2.7784504406519571e-05  3.0964622176844602e-03
-   16  2.0348857368272344e-03 -1.2415016403438035e-04  3.3675205082768264e-02
-   17 -2.7089567432503573e-03 -8.1788662401764114e-04  2.6578231998707000e-02
-   18  1.3589859834844708e-03 -4.1869284075205220e-04  3.2382362625908642e-02
-   19  3.0575998684987288e-04 -8.5615261685325662e-04  2.6300918365183709e-02
-   20  1.1482382726252934e-03 -2.6885864047674810e-04  3.4194853807722227e-02
-   21 -1.2680891310638792e-03  4.9877122346546224e-04  6.1154164862940522e-03
-   22  2.3554944751877239e-03  9.6068685589809240e-04  2.9850487598267021e-02
-   23  4.9959016556344255e-05  3.1673714227676154e-04  3.6614479955528112e-03
-   24  2.0192733259603224e-03  6.3534849028298586e-04  3.3638585635455770e-02
-   25 -2.3608563352120440e-03 -7.9888874314704573e-04  2.8465994834953420e-02
-   26  1.1497012824588045e-03 -3.9194609517719160e-04  3.2842627585052867e-02
-   27  2.9249325092285189e-04 -8.3077344407106490e-04  2.8308269755484061e-02
-   28  9.3417705968603434e-04 -2.3040781495046246e-04  3.4278191255128133e-02
-   29 -5.1652528928022433e-04 -1.4404353042683099e-04 -1.2462150875064196e-02
-   30  6.9288940488795056e-05 -1.4710788156542239e-04 -1.2875332550432134e-02
-   31  2.7644608312646314e-04 -1.4573040212804145e-04 -1.2258366148270420e-02
-   32  1.7312052380602741e-04 -1.3246895843787418e-04 -1.4250473428514749e-02
-   33 -4.2201273060055557e-04 -7.2075937985441126e-05 -1.3820153961237609e-02
-   34  3.3549241776445280e-05 -5.8924381188893540e-05 -1.3900492261623465e-02
-   35  2.6119063773180797e-04 -7.3975274656908229e-05 -1.3656086631368676e-02
-   36  1.2829798186810219e-04 -4.4146305074977061e-05 -1.4997087695938056e-02
-   37 -5.1375024459660446e-04  3.1547360096919003e-04 -1.2519087424517190e-02
-   38  6.8697062203298403e-05  2.9188411177715840e-04 -1.2927175328444108e-02
-   39  2.7540798890579723e-04  3.2141722701930868e-04 -1.2316658300695350e-02
-   40  1.7191884579607650e-04  2.4497453169432579e-04 -1.4294679585099376e-02
-   41 -3.8125358429239129e-04 -9.8073404256391242e-05 -1.4573806404289218e-02
-   42  2.5208712883704254e-05 -8.5394829999053783e-05 -1.4603558816163261e-02
-   43  2.4329125095209453e-04 -1.0051714757303402e-04 -1.4423743772053488e-02
-   44  1.1343880007511943e-04 -6.7742613833467660e-05 -1.5588639087563531e-02
+    1  2.4838446863956616e-04  1.9503788466163226e-03  3.0873204426106526e+00
+    2 -1.6652790386577014e-02  2.1026187422225435e-03  3.0749612676149027e+00
+    3  6.3785749905510489e-04  5.1125747155211397e-03  3.0896293958296610e+00
+    4 -1.6605063813845172e-02  6.0649203065202620e-03  3.0777057807809554e+00
+    5  4.2242722672783749e-02 -2.3050714764837649e-03 -1.5592361985876177e+00
+    6  5.5714273062644336e-02 -2.4127236407636803e-03 -1.5626351195306318e+00
+    7  4.2228746222139832e-02 -7.5938674438442760e-04 -1.5603926320308752e+00
+    8  5.5886401300220948e-02 -9.9610508505008281e-04 -1.5639533086735800e+00
+    9 -5.1076262902731158e-02 -2.3218371704875301e-03 -1.5701372726008263e+00
+   10 -3.3811561421032657e-02 -2.3397296264742822e-03 -1.5550660436064543e+00
+   11 -5.1350902848165500e-02 -5.0363947736732217e-04 -1.5712040115797643e+00
+   12 -3.4024564072498521e-02 -1.0255079009805182e-03 -1.5563752159595889e+00
+   13  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   14  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   15  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   16  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   17  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   18  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   19  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   20  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   21  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   22  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   23  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   24  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   25  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   26  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   27  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   28  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   29  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   30  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   31  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   32  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   33  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   34  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   35  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   36  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   37  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   38  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   39  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   40  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   41  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   42  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   43  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+   44  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 ...