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

Mix eta in Electrode package

Browse Source
This commit is contained in:
Ludwig Ahrens-Iwers
2023-11-13 14:41:13 +01:00
committed by Shern Tee
parent 857cc53923
commit a651697d2e
16 changed files with 303 additions and 78 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
doc/src/fix_electrode.rst
View File

@ -70,6 +70,8 @@ Syntax
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*
*eta* value = d_propname
d_propname = a custom double vector defined via fix property/atom
Examples
""""""""
@ -264,6 +266,11 @@ 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).
The keyword *eta* takes the name of a custom double vector defined via fix
property/atom. The values will be used instead of the standard eta value. The
property/atom fix must be for vector of double values and use the *ghost on*
option.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

34
examples/PACKAGES/electrode/madelung/data.eta Normal file
View File

@ -0,0 +1,34 @@
LAMMPS data file via write_data, version 24 Dec 2020, timestep = 0
4 atoms
3 atom types
0 1 xlo xhi
0 1 ylo yhi
-10 10 zlo zhi
Masses
1 196.966553
2 196.966553
3 1.0
Pair Coeffs # lj/cut/coul/long
1 0 0
2 0 0
3 0 0
Atoms # full
1 1 1 0.00 0.00 0.00 -2.00 # bottom electrode
2 2 2 0.00 0.00 0.00 2.00 # top electrode
3 3 3 0.50 0.00 0.00 -1.00 # bottom electrolyte
4 3 3 -0.50 0.00 0.00 1.00 # top electrolyte
ETA
1 2.0
2 2.0
3 0
4 0

34
examples/PACKAGES/electrode/madelung/data.eta_mix Normal file
View File

@ -0,0 +1,34 @@
LAMMPS data file via write_data, version 24 Dec 2020, timestep = 0
4 atoms
3 atom types
0 1 xlo xhi
0 1 ylo yhi
-10 10 zlo zhi
Masses
1 196.966553
2 196.966553
3 1.0
Pair Coeffs # lj/cut/coul/long
1 0 0
2 0 0
3 0 0
Atoms # full
1 1 1 0.00 0.00 0.00 -2.00 # bottom electrode
2 2 2 0.00 0.00 0.00 2.00 # top electrode
3 3 3 0.50 0.00 0.00 -1.00 # bottom electrolyte
4 3 3 -0.50 0.00 0.00 1.00 # top electrolyte
ETA
1 0.5
2 3.0
3 0
4 0

4
examples/PACKAGES/electrode/madelung/eval.py
View File

@ -1,7 +1,7 @@
#!/usr/env/python3
import sys
import os.path as op
import sys
def rel_error(out, ref):
@ -49,5 +49,5 @@ for label, ref, out in out_lines:
error = rel_error(out, ref)
lines.append(f"{label}: {out:.5f}, {error:.5f}\n")
with open("madelung.txt", 'a') as f:
with open("madelung.txt", "a") as f:
f.writelines(lines)

17
examples/PACKAGES/electrode/madelung/in.eta Normal file
View File

@ -0,0 +1,17 @@
atom_style full
units real
boundary p p f
kspace_style ewald/electrode 1.0e-8
kspace_modify slab 8.0 # ew3dc
pair_style lj/cut/coul/long 12
fix feta all property/atom d_eta ghost yes
read_data data.eta fix feta NULL ETA
include "settings_eta.mod"
fix conp bot electrode/conp 0 2 couple top 1 symm on eta d_eta write_inv inv.csv write_vec vec.csv
run 0

17
examples/PACKAGES/electrode/madelung/in.eta_cg Normal file
View File

@ -0,0 +1,17 @@
atom_style full
units real
boundary p p f
kspace_style ewald/electrode 1.0e-8
kspace_modify slab 8.0 # ew3dc
pair_style lj/cut/coul/long 12
fix feta all property/atom d_eta ghost yes
read_data data.eta_mix fix feta NULL ETA
include "settings_eta.mod"
fix conp bot electrode/conp 0 2 couple top 1 symm on algo cg 1e-6 eta d_eta
run 0

17
examples/PACKAGES/electrode/madelung/in.eta_mix Normal file
View File

@ -0,0 +1,17 @@
atom_style full
units real
boundary p p f
kspace_style ewald/electrode 1.0e-8
kspace_modify slab 8.0 # ew3dc
pair_style lj/cut/coul/long 12
fix feta all property/atom d_eta ghost on
read_data data.eta_mix fix feta NULL ETA
include "settings_eta.mod"
fix conp bot electrode/conp 0 2 couple top 1 symm on eta d_eta write_inv inv.csv write_vec vec.csv
run 0

83
examples/PACKAGES/electrode/madelung/plate_cap.py
View File

@ -3,7 +3,6 @@
import numpy as np
from scipy.special import erf
ETA = 2
SQRT2 = np.sqrt(2)
COULOMB = 332.06371 # Coulomb constant in Lammps 'real' units
QE2F = 23.060549
@ -17,14 +16,14 @@ def lattice(length):
return np.array(np.meshgrid(x, y)).T.reshape(-1, 2)
def a_element(r):
def a_element(r, eta):
"""Coulomb contribution of two Gaussians"""
return erf(ETA / SQRT2 * r) / r
return erf(eta * r) / r
def b_element(r, q):
def b_element(r, q, eta):
"""Coulomb contribution of a Gaussian with a point charge"""
return q * erf(ETA * r) / r
return q * erf(eta * r) / r
a = 1 # nearest neighbor distance i.e. lattice constant / sqrt(2)
@ -36,49 +35,54 @@ 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)
opposite_distances = np.sqrt(np.square(distances) + distance_plates**2)
# self interaction and within original box
A_11 = np.sqrt(2 / np.pi) * ETA
A_12 = erf(ETA * distance_plates / SQRT2) / distance_plates
for name, eta_elec in [("", [2.0, 2.0]), ("_eta_mix", [0.5, 3.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]
A_22 = np.sqrt(2 / np.pi) * eta_elec[1]
A_12 = erf(eta_mix * distance_plates) / distance_plates
# interaction with periodic images
A_11 += 4 * np.sum(a_element(distances))
A_12 += 4 * np.sum(a_element(opposite_distances))
A = np.array([[A_11, A_12], [A_12, A_11]])
inv = np.linalg.inv(A)
e = np.array([1, 1])
inv -= np.matmul(inv, np.matmul(np.outer(e, e), inv)) / np.dot(e, np.dot(inv, e))
# interaction with periodic images
A_11 += 4 * np.sum(a_element(distances, eta_elec[0] / SQRT2))
A_22 += 4 * np.sum(a_element(distances, eta_elec[1] / SQRT2))
A_12 += 4 * np.sum(a_element(opposite_distances, eta_mix))
A = np.array([[A_11, A_12], [A_12, A_22]])
inv = np.linalg.inv(A)
e = np.array([1, 1])
inv -= np.matmul(inv, np.matmul(np.outer(e, e), inv)) / np.dot(e, np.dot(inv, e))
# electrode-electrolyte interaction
b = []
for x in x_elec:
# electrode-electrolyte interaction
b = []
for x, eta in 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)
image_distances = np.sqrt(np.square(distances) + d ** 2)
bi += 4 * np.sum(b_element(image_distances, q))
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))
b.append(bi)
b = np.array(b)
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 = np.array([[elyt_11, elyt_12], [elyt_12, elyt_11]])
energy_elyt = 0.5 * np.dot(q_elyt, np.dot(elyt, q_elyt))
# 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 = np.array([[elyt_11, elyt_12], [elyt_12, elyt_11]])
energy_elyt = 0.5 * np.dot(q_elyt, np.dot(elyt, q_elyt))
# electrode charges and energy
q = np.dot(inv, v - b)
energy = COULOMB * (0.5 * np.dot(q, np.dot(A, q)) + np.dot(b, q) + energy_elyt)
# electrode charges and energy
q = np.dot(inv, v - b)
energy = COULOMB * (0.5 * np.dot(q, np.dot(A, q)) + np.dot(b, q) + energy_elyt)
print(
"length, energy / kcal/mol, q1 / e, q2 / e, inv11 / A, inv12 / A, b1 / e/A, b2 / e/A"
)
print(
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"
)
f.write(
", ".join(
[
str(LENGTH),
@ -91,4 +95,5 @@ print(
f"{b[1]:.8f}",
]
)
)
+ "\n"
)

19
examples/PACKAGES/electrode/madelung/settings_eta.mod Normal file
View File

@ -0,0 +1,19 @@
# distribute electrode atoms among all processors:
if "$(extract_setting(world_size) % 2) == 0" then "processors * * 2"
group bot type 1
group top type 2
# get electrode charges
variable q atom q
compute qbot bot reduce sum v_q
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"
thermo_style custom step pe c_qbot c_qtop

16
examples/PACKAGES/electrode/madelung/test.sh
View File

@ -7,17 +7,27 @@ if [ ! -f $lmpbin ]; then
fi
ref_out="plate_cap.csv"
if [ ! -f $ref_out ]; then
ref_mix_out="plate_cap_eta_mix.csv"
if [ ! -f $ref_out ] || [ ! -f $ref_mix_out ]; then
echo "Generating reference data"
python3 plate_cap.py > $ref_out
python3 plate_cap.py
fi
echo "Running Lammps inputs"
# w/o eta mixing
rm -rf madelung.txt && touch madelung.txt
for file in in.*; do
for file in in.eta in.ewald-ew3dc in.ewald-ew2d in.pppm-ew3dc in.cg; do
printf "\n$file\n" >> madelung.txt
rm -f out.csv inv.csv vec.csv
$lmpbin -i $file &> /dev/null
python3 eval.py $ref_out out.csv inv.csv vec.csv
done
# with eta mixing
for file in in.eta_mix in.eta_cg; do
printf "\n$file\n" >> madelung.txt
rm -f out.csv inv.csv vec.csv
$lmpbin -i $file &> /dev/null
python3 eval.py $ref_mix_out out.csv inv.csv vec.csv
done
cat madelung.txt

20
src/ELECTRODE/electrode_matrix.cpp
View File

@ -43,6 +43,7 @@ ElectrodeMatrix::ElectrodeMatrix(LAMMPS *lmp, int electrode_group, double eta) :
groupbit = group->bitmask[igroup];
ngroup = group->count(igroup);
this->eta = eta;
etaflag = false;
tfflag = false;
}
@ -72,6 +73,14 @@ void ElectrodeMatrix::setup_tf(const std::map<int, double> &tf_types)
/* ---------------------------------------------------------------------- */
void ElectrodeMatrix::setup_eta(int index)
{
etaflag = true;
eta_index = index;
}
/* ---------------------------------------------------------------------- */
void ElectrodeMatrix::compute_array(double **array, bool timer_flag)
{
// setting all entries of coulomb matrix to zero
@ -115,8 +124,6 @@ void ElectrodeMatrix::pair_contribution(double **array)
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double const etaij = eta * eta / sqrt(2.0 * eta * eta); // see mw ewald theory eq. (29)-(30)
// neighbor list will be ready because called from post_neighbor
inum = list->inum;
ilist = list->ilist;
@ -135,6 +142,7 @@ void ElectrodeMatrix::pair_contribution(double **array)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
double const eta_i = etaflag ? atom->dvector[eta_index][i] : eta;
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
@ -152,6 +160,9 @@ void ElectrodeMatrix::pair_contribution(double **array)
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
double const eta_j = etaflag ? atom->dvector[eta_index][j] : eta;
double const etaij = eta_i * eta_j / sqrt(eta_i * eta_i + eta_j * eta_j);
r = sqrt(rsq);
rinv = 1.0 / r;
aij = rinv;
@ -178,7 +189,10 @@ void ElectrodeMatrix::self_contribution(double **array)
const double preta = MY_SQRT2 / MY_PIS;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { array[mpos[i]][mpos[i]] += preta * eta - selfint; }
if (mask[i] & groupbit) {
double const eta_i = etaflag ? atom->dvector[eta_index][i] : eta;
array[mpos[i]][mpos[i]] += preta * eta_i - selfint;
}
}
/* ---------------------------------------------------------------------- */

3
src/ELECTRODE/electrode_matrix.h
View File

@ -30,6 +30,7 @@ class ElectrodeMatrix : protected Pointers {
ElectrodeMatrix(class LAMMPS *, int, double);
void setup(const std::unordered_map<tagint, int> &, class Pair *, class NeighList *);
void setup_tf(const std::map<int, double> &);
void setup_eta(int);
void compute_array(double **, bool);
int igroup;
@ -39,6 +40,8 @@ class ElectrodeMatrix : protected Pointers {
double **cutsq;
double g_ewald, eta;
bool tfflag;
bool etaflag;
int eta_index;
std::map<int, double> tf_types;
std::unordered_map<tagint, int> tag_to_iele;
bool assigned;

33
src/ELECTRODE/electrode_vector.cpp
View File

@ -29,6 +29,7 @@
#include "neigh_list.h"
#include "pair.h"
#include <cassert>
#include <cmath>
#include <exception>
@ -47,6 +48,7 @@ ElectrodeVector::ElectrodeVector(LAMMPS *lmp, int sensor_group, int source_group
source_grpbit = group->bitmask[source_group];
this->eta = eta;
tfflag = false;
etaflag = false;
kspace_time_total = 0;
pair_time_total = 0;
@ -93,6 +95,14 @@ void ElectrodeVector::setup_tf(const std::map<int, double> &tf_types)
/* ---------------------------------------------------------------------- */
void ElectrodeVector::setup_eta(int index)
{
etaflag = true;
eta_index = index;
}
/* ---------------------------------------------------------------------- */
void ElectrodeVector::compute_vector(double *vector)
{
MPI_Barrier(world);
@ -121,7 +131,6 @@ void ElectrodeVector::compute_vector(double *vector)
void ElectrodeVector::pair_contribution(double *vector)
{
double const etaij = eta * MY_ISQRT2;
double **x = atom->x;
double *q = atom->q;
int *type = atom->type;
@ -142,6 +151,7 @@ void ElectrodeVector::pair_contribution(double *vector)
double const xtmp = x[i][0];
double const ytmp = x[i][1];
double const ztmp = x[i][2];
double const eta_i = etaflag ? atom->dvector[eta_index][i] : eta;
int itype = type[i];
int *jlist = firstneigh[i];
int jnum = numneigh[i];
@ -158,18 +168,22 @@ void ElectrodeVector::pair_contribution(double *vector)
double const rsq = delx * delx + dely * dely + delz * delz;
int jtype = type[j];
if (rsq >= cutsq[itype][jtype]) continue;
double const eta_j = etaflag ? atom->dvector[eta_index][j] : eta;
double etaij;
if (i_in_sensor && j_in_sensor)
etaij = eta_i * eta_j / sqrt(eta_i * eta_i + eta_j * eta_j);
else if (i_in_sensor)
etaij = eta_i;
else {
assert(j_in_sensor);
etaij = eta_j;
}
double const r = sqrt(rsq);
double const rinv = 1.0 / r;
double aij = rinv;
aij *= ElectrodeMath::safe_erfc(g_ewald * r);
if (invert_source)
aij -= ElectrodeMath::safe_erfc(eta * r) * rinv;
else
aij -= ElectrodeMath::safe_erfc(etaij * r) * rinv;
if (i_in_sensor) {
vector[i] += aij * q[j];
//} else if (j_in_sensor) {
}
if (i_in_sensor) { vector[i] += aij * q[j]; }
if (j_in_sensor && (!invert_source || !i_in_sensor)) { vector[j] += aij * q[i]; }
}
}
@ -189,9 +203,10 @@ void ElectrodeVector::self_contribution(double *vector)
for (int ii = 0; ii < inum; ii++) {
int const i = ilist[ii];
double const eta_i = etaflag ? atom->dvector[eta_index][i] : eta;
bool const i_in_sensor = (mask[i] & groupbit);
bool const i_in_source = !!(mask[i] & source_grpbit) != invert_source;
if (i_in_sensor && i_in_source) vector[i] += (preta * eta - selfint) * q[i];
if (i_in_sensor && i_in_source) vector[i] += (preta * eta_i - selfint) * q[i];
}
}

3
src/ELECTRODE/electrode_vector.h
View File

@ -29,6 +29,7 @@ class ElectrodeVector : protected Pointers {
~ElectrodeVector() override;
void setup(class Pair *, class NeighList *, bool);
void setup_tf(const std::map<int, double> &);
void setup_eta(int);
void compute_vector(double *);
int igroup, source_group;
@ -39,6 +40,8 @@ class ElectrodeVector : protected Pointers {
double **cutsq;
double g_ewald, eta;
bool tfflag;
bool etaflag;
int eta_index;
std::map<int, double> tf_types;
class Pair *pair;
class NeighList *list;

36
src/ELECTRODE/fix_electrode_conp.cpp
View File

@ -97,6 +97,7 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) :
top_group = 0;
intelflag = false;
tfflag = false;
etaflag = false;
timer_flag = false;
update_time = 0;
@ -211,6 +212,19 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) :
qtotal = utils::numeric(FLERR, arg[iarg], false, lmp);
qtotal_var_style = VarStyle::CONST;
}
} else if ((strcmp(arg[iarg], "eta") == 0)) {
if (iarg + 2 > narg) error->all(FLERR, "Need two arguments after eta command");
etaflag = true;
int is_double, cols;
eta_index = atom->find_custom(arg[++iarg] + 2, is_double, cols);
if (eta_index == -1)
error->all(FLERR, "eta keyword requires name of previously defined property");
if (!is_double) error->all(FLERR, "eta keyword requires double-valued property/atom vector");
if (cols != 0) error->all(FLERR, "eta keyword requires property/atom vector not an array");
if (!atom->nextra_border)
error->all(FLERR,
"There is no fix with ghost on, but the eta keyword requires a property/atom "
"fix with ghost on");
// toggle parameters
} else if ((strcmp(arg[iarg], "etypes") == 0)) {
etypes_neighlists = utils::logical(FLERR, arg[++iarg], false, lmp);
@ -520,8 +534,10 @@ void FixElectrodeConp::setup_post_neighbor()
evscale = force->qe2f / force->qqrd2e;
elyt_vector->setup(pair, vec_neighlist, timer_flag);
if (etaflag) elyt_vector->setup_eta(eta_index);
if (need_elec_vector) {
elec_vector->setup(pair, mat_neighlist, timer_flag);
if (etaflag) elec_vector->setup_eta(eta_index);
if (tfflag) elec_vector->setup_tf(tf_types);
}
@ -556,7 +572,8 @@ void FixElectrodeConp::setup_post_neighbor()
if (etypes_neighlists) neighbor->build_one(mat_neighlist, 0);
auto array_compute = std::unique_ptr<ElectrodeMatrix>(new ElectrodeMatrix(lmp, igroup, eta));
array_compute->setup(tag_to_iele, pair, mat_neighlist);
if (tfflag) { array_compute->setup_tf(tf_types); }
if (etaflag) array_compute->setup_eta(eta_index);
if (tfflag) array_compute->setup_tf(tf_types);
array_compute->compute_array(elastance, timer_flag);
} // write_mat before proceeding
if (comm->me == 0 && write_mat) {
@ -1186,7 +1203,7 @@ double FixElectrodeConp::self_energy(int eflag)
// corrections to energy due to self interaction
double const qqrd2e = force->qqrd2e;
int const nlocal = atom->nlocal;
double const pre = eta / sqrt(MY_2PI) * qqrd2e;
double const pre = 1. / sqrt(MY_2PI) * qqrd2e;
int *mask = atom->mask;
int *type = atom->type;
double *q = atom->q;
@ -1194,7 +1211,8 @@ double FixElectrodeConp::self_energy(int eflag)
for (int i = 0; i < nlocal; i++) {
if (groupbit & mask[i]) {
double const q2 = q[i] * q[i];
double e = pre * q2;
double ieta = etaflag ? atom->dvector[eta_index][i] : eta;
double e = ieta * pre * q2;
if (tfflag && (groupbit & mask[i])) e += 0.5 * qqrd2e * q2 * tf_types[type[i]];
energy += e;
if (eflag) {
@ -1234,6 +1252,7 @@ double FixElectrodeConp::gausscorr(int eflag, bool fflag)
double xtmp = x[i][0];
double ytmp = x[i][1];
double ztmp = x[i][2];
double const eta_i = etaflag ? atom->dvector[eta_index][i] : eta;
int itype = type[i];
int *jlist = firstneigh[i];
int jnum = numneigh[i];
@ -1242,7 +1261,6 @@ double FixElectrodeConp::gausscorr(int eflag, bool fflag)
int const j = jlist[jj] & NEIGHMASK;
bool j_in_ele = groupbit & mask[j];
if (!(i_in_ele || j_in_ele)) continue;
double eta_ij = (i_in_ele && j_in_ele) ? eta / MY_SQRT2 : eta;
double delx = xtmp - x[j][0];
double dely = ytmp - x[j][1];
@ -1251,6 +1269,16 @@ double FixElectrodeConp::gausscorr(int eflag, bool fflag)
int jtype = type[j];
if (rsq < force->pair->cutsq[itype][jtype]) {
double const eta_j = etaflag ? atom->dvector[eta_index][j] : eta;
double eta_ij;
if (i_in_ele && j_in_ele)
eta_ij = eta_i * eta_j / sqrt(eta_i * eta_i + eta_j * eta_j);
else if (i_in_ele)
eta_ij = eta_i;
else {
assert(j_in_ele);
eta_ij = eta_j;
}
double r2inv = 1.0 / rsq;
double r = sqrt(rsq);
double erfc_etar = 0.;

2
src/ELECTRODE/fix_electrode_conp.h
View File

@ -138,6 +138,8 @@ class FixElectrodeConp : public Fix {
int get_top_group(); // used by ffield
int top_group; // used by ffield
bool tfflag;
int eta_index; // index of atom property for eta
bool etaflag; // eta specified as atom property
bool timer_flag;
std::map<int, double> tf_types;
// cg