From 49c84dbe1eb9f57c422f7d0e177df6031ba0d22e Mon Sep 17 00:00:00 2001 From: Ludwig Ahrens-Iwers Date: Mon, 1 Jul 2024 10:05:28 +0200 Subject: [PATCH] Bugfix virial with fix electrode --- .../PACKAGES/electrode/madelung/.gitignore | 1 + examples/PACKAGES/electrode/madelung/eval.py | 18 +++- .../PACKAGES/electrode/madelung/plate_cap.py | 98 +++++++++++++++-- .../PACKAGES/electrode/madelung/settings.mod | 6 +- src/ELECTRODE/fix_electrode_conp.cpp | 100 +++++++++++++++--- src/ELECTRODE/fix_electrode_conp.h | 3 +- 6 files changed, 197 insertions(+), 29 deletions(-) diff --git a/examples/PACKAGES/electrode/madelung/.gitignore b/examples/PACKAGES/electrode/madelung/.gitignore index 89d8d1a065..f222840fd6 100644 --- a/examples/PACKAGES/electrode/madelung/.gitignore +++ b/examples/PACKAGES/electrode/madelung/.gitignore @@ -1,2 +1,3 @@ *.csv *.txt +*.lammpstrj diff --git a/examples/PACKAGES/electrode/madelung/eval.py b/examples/PACKAGES/electrode/madelung/eval.py index feda0e384e..8f3675a741 100644 --- a/examples/PACKAGES/electrode/madelung/eval.py +++ b/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) diff --git a/examples/PACKAGES/electrode/madelung/plate_cap.py b/examples/PACKAGES/electrode/madelung/plate_cap.py index fcca166869..37aad47f51 100755 --- a/examples/PACKAGES/electrode/madelung/plate_cap.py +++ b/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): - 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)) + for j, (y, q) in enumerate(zip(x_elyt, q_elyt)): + bi += b_element(np.abs(y - x), q, eta) + bi += 4 * np.sum(b_element(image_distances[i][j], 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_12 = 1 / distance_elyt + 4 * np.sum(1 / image_elyt_distances[0][1]) 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") diff --git a/examples/PACKAGES/electrode/madelung/settings.mod b/examples/PACKAGES/electrode/madelung/settings.mod index aa1096ea81..bb5c8e42ae 100644 --- a/examples/PACKAGES/electrode/madelung/settings.mod +++ b/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 diff --git a/src/ELECTRODE/fix_electrode_conp.cpp b/src/ELECTRODE/fix_electrode_conp.cpp index a4b9d5b3b9..aaa7f2fb99 100644 --- a/src/ELECTRODE/fix_electrode_conp.cpp +++ b/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; } } - - double ecoul = 0.; - if (eflag) ecoul = -prefactor * erfc_etar; - - if (evflag) { - force->pair->ev_tally(i, j, nlocal, newton_pair, 0., ecoul, fpair, delx, dely, delz); + if (eflag) { + double ecoul = -prefactor * erfc_etar; + force->pair->ev_tally(i, j, nlocal, newton_pair, 0., ecoul, 0., 0., 0., 0.); } + 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]; + } + } + } +} diff --git a/src/ELECTRODE/fix_electrode_conp.h b/src/ELECTRODE/fix_electrode_conp.h index a1d7530bd1..b0b4a1fd46 100644 --- a/src/ELECTRODE/fix_electrode_conp.h +++ b/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 &, const std::vector> &); void read_from_file(const std::string &input_file, double **, const std::string &); void compute_sd_vectors();