diff --git a/doc/src/fix_electrode.rst b/doc/src/fix_electrode.rst index 3d543f08d2..5685482a71 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* .. parsed-literal:: @@ -68,6 +68,8 @@ 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* Examples """""""" @@ -249,6 +251,19 @@ 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). + Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/examples/PACKAGES/electrode/graph-il/in.nonneut b/examples/PACKAGES/electrode/graph-il/in.nonneut new file mode 100644 index 0000000000..3027faf0cc --- /dev/null +++ b/examples/PACKAGES/electrode/graph-il/in.nonneut @@ -0,0 +1,31 @@ +# electrodes, overall non-neutral +# for graphene-ionic liquid supercapacitor + +boundary p p f # slab calculation +include settings.mod # styles, groups, computes and fixes +kspace_modify slab 3.0 + +# make an extra anion +variable qmodify index 1 +variable vdiff index 1 +# these values can be changed, e.g. ${LAMMPS_EXECUTABLE} -i in.nonneut -var qmodify 0.9 -var vdiff -2 +# fix electrode/* needs equal style variables: +variable qmodify_equal equal v_qmodify +variable vdiff_equal equal v_vdiff + +create_atoms 4 single 16.1 17.2 0 +set atom 3777 charge $(-v_qmodify) + +fix c top electrode/conp v_vdiff_equal 1.979 couple bot 0 etypes on qtotal v_qmodify_equal +# to test electrode/thermo: +# variable vbot equal 0 +# fix c top electrode/thermo v_vdiff_equal 1.979 couple bot v_vbot etypes on qtotal v_qmodify_equal temp 310 100 12309 # symm on + +variable dv equal f_c[1]-f_c[2] + +variable qelec equal c_qbot+c_qtop +compute qall all reduce sum v_q # total system charge + +thermo 50 +thermo_style custom step temp c_ctemp epair etotal c_qbot c_qtop f_c[1] f_c[2] v_dv v_qelec c_qall +run 500 diff --git a/examples/PACKAGES/electrode/graph-il/in.nonneut2 b/examples/PACKAGES/electrode/graph-il/in.nonneut2 new file mode 100644 index 0000000000..c73eee5db3 --- /dev/null +++ b/examples/PACKAGES/electrode/graph-il/in.nonneut2 @@ -0,0 +1,33 @@ +# electrodes, overall non-neutral +# for graphene-ionic liquid supercapacitor + +boundary p p f # slab calculation +include settings.mod # styles, groups, computes and fixes +kspace_modify slab 3.0 + +# make an extra anion +variable qmodify index 1 +variable vdiff index 1 +# these values can be changed, e.g. ${LAMMPS_EXECUTABLE} -i in.nonneut -var qmodify 0.9 -var vdiff -2 +create_atoms 4 single 16.1 17.2 0 +set atom 3777 charge $(-v_qmodify) + +fix c top electrode/conp v_vtop 1.979 couple bot v_vbot etypes on # symm on +# to test electrode/thermo: +#fix c top electrode/thermo v_vtop 1.979 couple bot v_vbot etypes on temp 310 100 12309 symm off +#but it looks like update_psi() is hardwired to impose electroneutrality + +variable csum equal f_c[1][2]+f_c[1][3]+f_c[2][2]+f_c[2][3] +variable cdiff equal 0.5*(f_c[1][2]-f_c[2][3]) +variable qdeficit equal v_qmodify-f_c[1][1]-f_c[2][1]-v_cdiff*v_vdiff +variable vphi equal v_qdeficit/v_csum +variable vtop equal v_vphi+0.5*v_vdiff +variable vbot equal v_vphi-0.5*v_vdiff +variable dv equal v_vtop-v_vbot + +variable qelec equal c_qbot+c_qtop +compute qall all reduce sum v_q # total system charge + +thermo 50 +thermo_style custom step temp c_ctemp epair etotal c_qbot c_qtop f_c[1] f_c[2] v_dv v_qelec c_qall +run 500 diff --git a/src/ELECTRODE/fix_electrode_conp.cpp b/src/ELECTRODE/fix_electrode_conp.cpp index 33632f9099..a7e6cbe922 100644 --- a/src/ELECTRODE/fix_electrode_conp.cpp +++ b/src/ELECTRODE/fix_electrode_conp.cpp @@ -103,20 +103,23 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) : 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(1, igroup); group_bits = std::vector(1, groupbit); group_psi_var_names = std::vector(1); group_psi_var_styles = std::vector(1, VarStyle::CONST); - group_psi = std::vector(1); + group_psi_const = std::vector(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; @@ -132,12 +135,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, fmt::format("Algorithm can be set once, only")); @@ -195,6 +198,19 @@ 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); + } 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; + } // toggle parameters } else if ((strcmp(arg[iarg], "etypes") == 0)) { etypes_neighlists = utils::logical(FLERR, arg[++iarg], false, lmp); @@ -208,6 +224,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(groups.size()); // union of all coupled groups std::string union_group = "conp_group"; std::string group_cmd = union_group + " union"; @@ -225,6 +247,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); @@ -396,7 +419,8 @@ void FixElectrodeConp::init() for (char *fix_id : integrate_ids) error->warning( FLERR, - "Electrode atoms are integrated by fix {}, but fix electrode is using a matrix method. For " + "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); } @@ -475,6 +499,7 @@ void FixElectrodeConp::setup_post_neighbor() // get equal-style variable ids: group_psi_var_ids = std::vector(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); @@ -485,6 +510,16 @@ void FixElectrodeConp::setup_post_neighbor() fmt::format("Variable '{}' for fix electrode is not equal-style", var_name)); 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, fmt::format("Variable '{}' for fix electrode does not exist", var_name)); + if (!input->variable->equalstyle(var_id)) + error->all(FLERR, + fmt::format("Variable '{}' for fix electrode is not equal-style", var_name)); + qtotal_var_id = var_id; + } // pair and list setups: @@ -834,6 +869,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); @@ -936,12 +973,22 @@ std::vector FixElectrodeConp::gather_ngroup(std::vector 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 FixElectrodeConp::constraint_correction(std::vector x) { - return constraint_projection(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; } /* ---------------------------------------------------------------------- @@ -950,7 +997,7 @@ std::vector FixElectrodeConp::constraint_correction(std::vector std::vector FixElectrodeConp::constraint_projection(std::vector 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); @@ -1008,13 +1055,28 @@ std::vector FixElectrodeConp::times_elastance(std::vector 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); @@ -1070,6 +1132,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]; } /* ---------------------------------------------------------------------- */ diff --git a/src/ELECTRODE/fix_electrode_conp.h b/src/ELECTRODE/fix_electrode_conp.h index 1289d96281..689a44053c 100644 --- a/src/ELECTRODE/fix_electrode_conp.h +++ b/src/ELECTRODE/fix_electrode_conp.h @@ -71,7 +71,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, std::vector){}; @@ -80,6 +80,7 @@ class FixElectrodeConp : public Fix { virtual void compute_macro_matrices(); std::vector ele_ele_interaction(const std::vector &); std::vector group_psi; + std::vector group_psi_const; // needed to undo qtotal psi updates std::vector group_bits; std::vector groups; int num_of_groups; @@ -101,6 +102,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; @@ -143,6 +148,9 @@ class FixElectrodeConp : public Fix { std::vector gather_ngroup(std::vector); std::vector gather_elevec_local(ElectrodeVector *); void set_charges(std::vector); + // qtotal + double macro_capacitance_sum; + void update_psi_qtotal(); // fix-specific electrode ID storage system: diff --git a/src/ELECTRODE/fix_electrode_thermo.cpp b/src/ELECTRODE/fix_electrode_thermo.cpp index 52c0a3ce4c..3883aa40a2 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]); }