ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

enable symm for conq and thermo

Browse Source
This commit is contained in:
Shern Ren Tee
2022-11-22 14:39:25 +10:00
committed by Shern Tee
parent 2efd485bd0
commit 32a3fc21bc
4 changed files with 33 additions and 29 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
examples/PACKAGES/electrode/graph-il/in.conq
View File

@ -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

2
examples/PACKAGES/electrode/graph-il/in.thermo
View File

@ -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

20
src/ELECTRODE/fix_electrode_conq.cpp
View File

@ -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<double> 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<double> 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<double> b, std::vector<do
int const n = b.size();
auto a = ele_ele_interaction(q_local);
auto psi_sums = std::vector<double>(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]);
}

34
src/ELECTRODE/fix_electrode_thermo.cpp
View File

@ -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];
}