diff --git a/examples/PACKAGES/electrode/graph-il/in.conq b/examples/PACKAGES/electrode/graph-il/in.conq index 9706ba3b9b..7f315a4b12 100644 --- a/examples/PACKAGES/electrode/graph-il/in.conq +++ b/examples/PACKAGES/electrode/graph-il/in.conq @@ -5,8 +5,10 @@ boundary p p f # slab calculation include settings.mod # styles, groups, computes and fixes kspace_modify slab 3.0 -fix conq bot electrode/conq -1.0 1.979 couple top 1.0 etypes on # conq doesn't take symm option +fix conq bot electrode/conq -1.0 1.979 couple top 1.0 etypes on # symm on +variable dv equal f_conq[2]-f_conq[1] +# symm on and off give different electrode potentials, but identical potential difference thermo 50 -thermo_style custom step temp c_ctemp epair etotal c_qbot c_qtop f_conq[1] f_conq[2] +thermo_style custom step temp c_ctemp epair etotal c_qbot c_qtop f_conq[1] f_conq[2] v_dv run 500 diff --git a/examples/PACKAGES/electrode/graph-il/in.thermo b/examples/PACKAGES/electrode/graph-il/in.thermo index 253dafc138..eed574bff6 100644 --- a/examples/PACKAGES/electrode/graph-il/in.thermo +++ b/examples/PACKAGES/electrode/graph-il/in.thermo @@ -6,7 +6,7 @@ include settings.mod # styles, groups, computes and fixes kspace_modify slab 3.0 unfix nvt # remove NVT thermostat included from "settings.mod" -fix conpthermo bot electrode/thermo -1.0 1.979 couple top 1.0 etypes on temp 500 100 7 # temp tau rng +fix conpthermo bot electrode/thermo -1.0 1.979 couple top 1.0 etypes on temp 500 100 7 symm on # temp tau rng # to compare to regular constant potential, switch previous line to this: #fix conp bot electrode/conp -1.0 1.979 couple top 1.0 etypes on symm on fix nve electrolyte nve diff --git a/src/ELECTRODE/fix_electrode_conq.cpp b/src/ELECTRODE/fix_electrode_conq.cpp index bad1500edf..ce2ba605f7 100644 --- a/src/ELECTRODE/fix_electrode_conq.cpp +++ b/src/ELECTRODE/fix_electrode_conq.cpp @@ -31,19 +31,25 @@ FixElectrodeConq::FixElectrodeConq(LAMMPS *lmp, int narg, char **arg) : // copy const-style values across because update_psi will change group_psi group_q = group_psi; - if (symm) error->all(FLERR, "Keyword symm on not allowed in electrode/conq"); + if (symm) { + if (num_of_groups == 1) + error->all(FLERR, "Keyword symm on not allowed in electrode/conq with only one electrode"); + error->warning(FLERR, + "Fix electrode/conq with keyword symm ignores the charge setting for the last " + "electrode listed"); + } } void FixElectrodeConq::update_psi() { - for (int g = 0; g < num_of_groups; g++) { + int const numsymm = num_of_groups - ((symm) ? 1 : 0); + for (int g = 0; g < numsymm; g++) { if (group_psi_var_styles[g] == VarStyle::CONST) continue; group_q[g] = input->variable->compute_equal(group_psi_var_ids[g]); } if (algo == Algo::MATRIX_INV) { - std::vector group_remainder_q(num_of_groups); - for (int g = 0; g < num_of_groups; g++) { group_remainder_q[g] = group_q[g] - sb_charges[g]; } - + std::vector group_remainder_q(num_of_groups, 0.); + for (int g = 0; g < numsymm; g++) { group_remainder_q[g] = group_q[g] - sb_charges[g]; } for (int g = 0; g < num_of_groups; g++) { double vtmp = 0; for (int h = 0; h < num_of_groups; h++) { @@ -98,9 +104,7 @@ void FixElectrodeConq::recompute_potential(std::vector b, std::vector(num_of_groups, 0); - for (int i = 0; i < n; i++) { - psi_sums[iele_to_group_local[i]] += (a[i] + b[i]) / evscale; - } + for (int i = 0; i < n; i++) { psi_sums[iele_to_group_local[i]] += (a[i] + b[i]) / evscale; } MPI_Allreduce(MPI_IN_PLACE, &psi_sums.front(), num_of_groups, MPI_DOUBLE, MPI_SUM, world); for (int g = 0; g < num_of_groups; g++) group_psi[g] = psi_sums[g] / group->count(groups[g]); } diff --git a/src/ELECTRODE/fix_electrode_thermo.cpp b/src/ELECTRODE/fix_electrode_thermo.cpp index a1ad9897bd..3bdf446f53 100644 --- a/src/ELECTRODE/fix_electrode_thermo.cpp +++ b/src/ELECTRODE/fix_electrode_thermo.cpp @@ -42,7 +42,6 @@ FixElectrodeThermo::FixElectrodeThermo(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR, "Number of electrodes != two in electrode/thermo"); if (group_psi_var_styles[0] != group_psi_var_styles[1]) error->all(FLERR, "Potentials in electrode/thermo must have same style"); - if (symm) error->all(FLERR, "Keyword symm on not allowed in electrode/thermo"); if (algo != Algo::MATRIX_INV) error->all(FLERR, "Algorithm not allowed in electrode/thermo"); if (thermo_time < SMALL) error->all(FLERR, "Keyword temp not set or zero in electrode/thermo"); @@ -62,9 +61,12 @@ FixElectrodeThermo::~FixElectrodeThermo() void FixElectrodeThermo::compute_macro_matrices() { FixElectrodeConp::compute_macro_matrices(); - vac_cap = (macro_capacitance[0][0] * macro_capacitance[1][1] - - macro_capacitance[0][1] * macro_capacitance[0][1]) / - (macro_capacitance[0][0] + macro_capacitance[1][1] + 2 * macro_capacitance[0][1]); + if (symm) + vac_cap = macro_capacitance[0][0]; + else + vac_cap = (macro_capacitance[0][0] * macro_capacitance[1][1] - + macro_capacitance[0][1] * macro_capacitance[0][1]) / + (macro_capacitance[0][0] + macro_capacitance[1][1] + 2 * macro_capacitance[0][1]); } /* ----------------------------------------------------------------------- */ @@ -91,14 +93,11 @@ void FixElectrodeThermo::update_psi() double const dt = update->dt; // group_q_eff is charge that corresponds to potential after previous step - double group_q_eff[NUM_GROUPS] = {0., 0.}; - for (int g = 0; g < NUM_GROUPS; g++) { group_q_eff[g] = group_q_old[g] - sb_charges[g]; } - double group_psi_old[NUM_GROUPS] = {0., 0.}; - for (int g = 0; g < NUM_GROUPS; g++) { - double vtmp = 0; - for (int h = 0; h < NUM_GROUPS; h++) { vtmp += macro_elastance[g][h] * group_q_eff[h]; } - group_psi_old[g] = vtmp; - } + double const group_q_eff[NUM_GROUPS] = {group_q_old[0] - sb_charges[0], + (symm) ? 0. : group_q_old[1] - sb_charges[1]}; + double const group_psi_old[NUM_GROUPS] = { + macro_elastance[0][0] * group_q_eff[0] + macro_elastance[0][1] * group_q_eff[1], + macro_elastance[1][0] * group_q_eff[0] + macro_elastance[1][1] * group_q_eff[1]}; double const delta_psi = group_psi_old[1] - group_psi_old[0]; // target potential difference from input parameters @@ -113,11 +112,10 @@ void FixElectrodeThermo::update_psi() thermo_random->gaussian(); double const group_remainder_q[NUM_GROUPS] = {-delta_charge - sb_charges[0], - delta_charge - sb_charges[1]}; + (symm) ? 0. : delta_charge - sb_charges[1]}; - for (int g = 0; g < NUM_GROUPS; g++) { - double vtmp = 0; - for (int h = 0; h < NUM_GROUPS; h++) { vtmp += macro_elastance[g][h] * group_remainder_q[h]; } - group_psi[g] = vtmp; - } + group_psi[0] = + macro_elastance[0][0] * group_remainder_q[0] + macro_elastance[0][1] * group_remainder_q[1]; + group_psi[1] = + macro_elastance[1][0] * group_remainder_q[0] + macro_elastance[1][1] * group_remainder_q[1]; }