Bugfix virial with fix electrode

examples/PACKAGES/electrode/madelung/.gitignore

@@ -1,2 +1,3 @@
 *.csv
 *.txt
+*.lammpstrj

examples/PACKAGES/electrode/madelung/eval.py

@@ -17,14 +17,22 @@ q_ref = float(ref_line[3])
 inv11_ref = float(ref_line[4])
 inv12_ref = float(ref_line[5])
 b1_ref = float(ref_line[6])
+felec1_ref = float(ref_line[8])
+felyt1_ref = float(ref_line[10])
+press_ref = float(ref_line[12])
 
 # out.csv
 with open(sys.argv[2]) as f:
     out_line = f.readlines()[-1].split(", ")
 e_out = float(out_line[0])
 q_out = float(out_line[1])
+press_out = float(out_line[2])
 
-out_lines = [("energy", e_ref, e_out), ("charge", q_ref, q_out)]
+out_lines = [
+    ("energy", e_ref, e_out),
+    ("charge", q_ref, q_out),
+    ("pressure", press_ref, press_out),
+]
 
 # vec.csv
 vec_file = "vec.csv"
@@ -44,6 +52,14 @@ if op.isfile(inv_file):
     inv12_out = float(inv_line[1])
     out_lines.append(("inv11", inv11_ref, inv11_out))
 
+# forces.lammpstrj
+force_file = "forces.lammpstrj"
+with open(force_file) as f:
+    lines = f.readlines()[9:]
+    for name, i, f_ref in [("felec1", "1", felec1_ref), ("felyt1", "3", felyt1_ref)]:
+        f_out = next(float(y[3]) for x in lines if (y := x.split()) and y[0] == i)
+        out_lines.append((name, f_ref, f_out))
+
 lines = []
 for label, ref, out in out_lines:
     error = rel_error(out, ref)

examples/PACKAGES/electrode/madelung/plate_cap.py

@@ -1,12 +1,17 @@
 #!/usr/bin/env python3
 
+import time
+
 import numpy as np
 from scipy.special import erf
 
 SQRT2 = np.sqrt(2)
+SQRTPI_INV = 1 / np.sqrt(np.pi)
 COULOMB = 332.06371  #  Coulomb constant in Lammps 'real' units
 QE2F = 23.060549
+NKTV2P = 68568.415  # pressure in 'real' units
 LENGTH = 10000  # convergence parameter
+LZ = 20
 
 
 def lattice(length):
@@ -26,6 +31,25 @@ def b_element(r, q, eta):
     return q * erf(eta * r) / r
 
 
+def force_gauss(r, qq, eta):
+    etar = eta * r
+    return (qq / np.square(r)) * (
+        erf(etar) - 2 * etar * SQRTPI_INV * np.exp(-np.square(etar))
+    )
+
+
+def force_point(r, qq):
+    return qq / np.square(r)
+
+
+def force_component(dx, d, qq, eta=None):
+    if eta:
+        return np.sum(dx / d * force_gauss(d, qq, eta))
+    else:
+        return np.sum(dx / d * force_point(d, qq))
+
+
+time_start = time.perf_counter()
 a = 1  # nearest neighbor distance i.e. lattice constant / sqrt(2)
 x_elec = [-2, 2]
 x_elyt = [-1, 1]
@@ -36,8 +60,20 @@ 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)
+image_distances = []
+for x in x_elec:
+    image_distances.append([])
+    for y in x_elyt:
+        image_distances[-1].append(np.sqrt(np.square(distances) + np.abs(y - x) ** 2))
+image_elyt_distances = [[None for _ in range(len(x_elyt))] for _ in range(len(x_elyt))]
+for i, (xi, qi) in enumerate(zip(x_elyt, q_elyt)):
+    for j, (xj, qj) in list(enumerate(zip(x_elyt, q_elyt)))[i + 1 :]:
+        image_elyt_distances[i][j] = np.sqrt(
+            np.square(distances) + np.abs(xj - xi) ** 2
+        )
 
 for name, eta_elec in [("", [2.0, 2.0]), ("_eta_mix", [0.5, 3.0])]:
+    # for name, eta_elec in [("", [2.0, 2.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]
@@ -55,22 +91,18 @@ for name, eta_elec in [("", [2.0, 2.0]), ("_eta_mix", [0.5, 3.0])]:
 
     # electrode-electrolyte interaction
     b = []
-    for x, eta in zip(x_elec, eta_elec):
+    for i, (x, eta) in enumerate(zip(x_elec, eta_elec)):
         bi = 0
-        for y, q in zip(x_elyt, q_elyt):
+        for j, (y, q) in enumerate(zip(x_elyt, q_elyt)):
-            d = abs(y - x)
+            bi += b_element(np.abs(y - x), q, eta)
-            bi += b_element(d, q, eta)
+            bi += 4 * np.sum(b_element(image_distances[i][j], 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(
+    elyt_12 = 1 / distance_elyt + 4 * np.sum(1 / image_elyt_distances[0][1])
-        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))
 
@@ -78,9 +110,48 @@ for name, eta_elec in [("", [2.0, 2.0]), ("_eta_mix", [0.5, 3.0])]:
     q = np.dot(inv, v - b)
     energy = COULOMB * (0.5 * np.dot(q, np.dot(A, q)) + np.dot(b, q) + energy_elyt)
 
+    # forces in out-of-plane direction
+    f_elec = np.zeros(len(x_elec))
+    f_elyt = np.zeros(len(x_elyt))
+    # electrode-electrode
+    dx = x_elec[1] - x_elec[0]
+    fij_box = force_component(dx, np.abs(dx), q[0] * q[1], eta_mix)
+    fij_img = 4 * force_component(dx, opposite_distances, q[0] * q[1], eta_mix)
+    f_elec[0] -= fij_box + fij_img
+    f_elec[1] += fij_box + fij_img
+    # electrode-electrolyte
+    for i, (xi, qi, etai) in enumerate(zip(x_elec, q, eta_elec)):
+        for j, (xj, qj) in enumerate(zip(x_elyt, q_elyt)):
+            dx = xj - xi
+            fij_box = force_component(dx, np.abs(dx), qi * qj, etai)
+            fij_img = 4 * force_component(dx, image_distances[i][j], qi * qj, etai)
+            f_elec[i] -= fij_box + fij_img
+            f_elyt[j] += fij_box + fij_img
+    # electrolyte-electrolyte
+    for i, (xi, qi) in enumerate(zip(x_elyt, q_elyt)):
+        for j, (xj, qj) in list(enumerate(zip(x_elyt, q_elyt)))[i + 1 :]:
+            dx = xj - xi
+            fij_box = force_component(dx, np.abs(dx), qi * qj)
+            fij_img = 4 * force_component(dx, image_elyt_distances[i][j], qi * qj)
+            f_elyt[i] -= fij_img + fij_box
+            f_elyt[j] += fij_img + fij_box
+    # force units
+    assert np.abs(np.sum(f_elec) + np.sum(f_elyt)) < 1e-8
+    f_elec *= COULOMB
+    f_elyt *= COULOMB
+
+    # Virial
+    volume = a**2 * LZ
+    virial = 0.0
+    for x, f in [(x_elec, f_elec), (x_elyt, f_elyt)]:
+        virial += np.dot(x, f)
+    pressure = NKTV2P * virial / volume
+
     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"
+            "length, energy / kcal/mol, q1 / e, q2 / e, inv11 / A, inv12 / A"
+            + ", b1 / e/A, b2 / e/A, felec1 / kcal/mol/A, felec2 / kcal/mol/A"
+            + ", felyt1 / kcal/mol/A, felyt2 / kcal/mol/A, press\n"
         )
         f.write(
             ", ".join(
@@ -93,7 +164,14 @@ for name, eta_elec in [("", [2.0, 2.0]), ("_eta_mix", [0.5, 3.0])]:
                     f"{inv[0, 1]:.10f}",
                     f"{b[0]:.8f}",
                     f"{b[1]:.8f}",
+                    f"{f_elec[0]:.5f}",
+                    f"{f_elec[1]:.5f}",
+                    f"{f_elyt[0]:.5f}",
+                    f"{f_elyt[1]:.5f}",
+                    f"{pressure:.2f}",
                 ]
             )
             + "\n"
         )
+time_end = time.perf_counter()
+print(f"{time_end - time_start:0.4f} seconds")

examples/PACKAGES/electrode/madelung/settings.mod

@@ -19,4 +19,8 @@ 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"
+compute press all pressure NULL virial
+variable p3 equal c_press[3]
+fix fxprint all print 1 "${vpe}, ${charge}, ${p3}" file "out.csv"
+
+dump dump_forces all custom 1 forces.lammpstrj id fx fy fz

src/ELECTRODE/fix_electrode_conp.cpp

@@ -89,6 +89,9 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) :
   extvector = 0;
   extarray = 0;
 
+  virial_global_flag = 1;    // use virials of this fix
+  thermo_virial = 1;         // set vflags for v_tally
+
   bool default_algo = true;
   algo = Algo::MATRIX_INV;
   matrix_algo = true;
@@ -642,13 +645,15 @@ void FixElectrodeConp::setup_post_neighbor()
 
 /* ---------------------------------------------------------------------- */
 
-void FixElectrodeConp::setup_pre_reverse(int eflag, int /*vflag*/)
+void FixElectrodeConp::setup_pre_reverse(int eflag, int vflag)
 {
+  if (pair->did_tally_callback())
+    error->warning(FLERR,
+                   "Computation of Virials in fix electrode is incompatible with TALLY package");
   // correct forces for initial timestep
-  gausscorr(eflag, true);
+  ev_init(eflag, vflag);
+  gausscorr(eflag, vflag, true);
   self_energy(eflag);
-  // potential_energy(eflag); // not always part of the energy, depending on ensemble, therefore
-  // removed
 }
 
 /* ---------------------------------------------------------------------- */
@@ -775,12 +780,11 @@ void FixElectrodeConp::pre_force(int)
 
 /* ---------------------------------------------------------------------- */
 
-void FixElectrodeConp::pre_reverse(int eflag, int /*vflag*/)
+void FixElectrodeConp::pre_reverse(int eflag, int vflag)
 {
-  gausscorr(eflag, true);
+  ev_init(eflag, vflag);
+  gausscorr(eflag, vflag, true);
   self_energy(eflag);
-  //potential_energy(eflag); // not always part of the energy, depending on ensemble, therefore
-  // removed
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1228,11 +1232,10 @@ double FixElectrodeConp::self_energy(int eflag)
 
 /* ---------------------------------------------------------------------- */
 
-double FixElectrodeConp::gausscorr(int eflag, bool fflag)
+double FixElectrodeConp::gausscorr(int eflag, int vflag, bool fflag)
 {
   // correction to short range interaction due to eta
 
-  int evflag = pair->evflag;
   double const qqrd2e = force->qqrd2e;
   int const nlocal = atom->nlocal;
   int *mask = atom->mask;
@@ -1298,13 +1301,11 @@ double FixElectrodeConp::gausscorr(int eflag, bool fflag)
             f[j][2] -= delz * fpair;
           }
         }
-
+        if (eflag) {
-        double ecoul = 0.;
+          double ecoul = -prefactor * erfc_etar;
-        if (eflag) ecoul = -prefactor * erfc_etar;
+          force->pair->ev_tally(i, j, nlocal, newton_pair, 0., ecoul, 0., 0., 0., 0.);
-
-        if (evflag) {
-          force->pair->ev_tally(i, j, nlocal, newton_pair, 0., ecoul, fpair, delx, dely, delz);
         }
+        if (vflag) v_tally(i, j, nlocal, newton_pair, fpair, delx, dely, delz);
       }
     }
   }
@@ -1625,3 +1626,70 @@ void FixElectrodeConp::unpack_forward_comm(int n, int first, double *buf)
   int const last = first + n;
   for (int i = first, m = 0; i < last; i++) atom->q[i] = buf[m++];
 }
+
+/* ----------------------------------------------------------------------
+   Tally virial of pair interactions in pre_reverse. This cannot be done with pair->ev_tally()
+   because compute_fdotr is called before pre_reverse, i.e. Virials need to be tallied even if fdotr
+   is used.
+------------------------------------------------------------------------- */
+
+void FixElectrodeConp::v_tally(int i, int j, int nlocal, int newton_pair, double fpair, double delx,
+                               double dely, double delz)
+{
+  double v[6];
+  if (vflag_either) {
+    v[0] = delx * delx * fpair;
+    v[1] = dely * dely * fpair;
+    v[2] = delz * delz * fpair;
+    v[3] = delx * dely * fpair;
+    v[4] = delx * delz * fpair;
+    v[5] = dely * delz * fpair;
+
+    if (vflag_global) {
+      if (newton_pair) {
+        virial[0] += v[0];
+        virial[1] += v[1];
+        virial[2] += v[2];
+        virial[3] += v[3];
+        virial[4] += v[4];
+        virial[5] += v[5];
+      } else {
+        if (i < nlocal) {
+          virial[0] += 0.5 * v[0];
+          virial[1] += 0.5 * v[1];
+          virial[2] += 0.5 * v[2];
+          virial[3] += 0.5 * v[3];
+          virial[4] += 0.5 * v[4];
+          virial[5] += 0.5 * v[5];
+        }
+        if (j < nlocal) {
+          virial[0] += 0.5 * v[0];
+          virial[1] += 0.5 * v[1];
+          virial[2] += 0.5 * v[2];
+          virial[3] += 0.5 * v[3];
+          virial[4] += 0.5 * v[4];
+          virial[5] += 0.5 * v[5];
+        }
+      }
+    }
+
+    if (vflag_atom) {
+      if (newton_pair || i < nlocal) {
+        vatom[i][0] += 0.5 * v[0];
+        vatom[i][1] += 0.5 * v[1];
+        vatom[i][2] += 0.5 * v[2];
+        vatom[i][3] += 0.5 * v[3];
+        vatom[i][4] += 0.5 * v[4];
+        vatom[i][5] += 0.5 * v[5];
+      }
+      if (newton_pair || j < nlocal) {
+        vatom[j][0] += 0.5 * v[0];
+        vatom[j][1] += 0.5 * v[1];
+        vatom[j][2] += 0.5 * v[2];
+        vatom[j][3] += 0.5 * v[3];
+        vatom[j][4] += 0.5 * v[4];
+        vatom[j][5] += 0.5 * v[5];
+      }
+    }
+  }
+}

src/ELECTRODE/fix_electrode_conp.h

@@ -119,10 +119,11 @@ class FixElectrodeConp : public Fix {
   void create_taglist();
   void invert();
   void symmetrize();
-  double gausscorr(int, bool);
+  double gausscorr(int, int, bool);
   void update_charges();
   double potential_energy();
   double self_energy(int);
+  void v_tally(int, int, int, int, double, double, double, double);
   void write_to_file(FILE *, const std::vector<tagint> &, const std::vector<std::vector<double>> &);
   void read_from_file(const std::string &input_file, double **, const std::string &);
   void compute_sd_vectors();