Skip to content

Efield lepton updates #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
gsalkuin merged 4 commits into from
Jan 15, 2025
Merged

Efield lepton updates #3

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
7b9f7be
replace vecs with arrays since size known at compile-time; make phi a…
srtee Jan 14, 2025
1f75330
deduplicate force and torque calculations
srtee Jan 14, 2025
276b8d9
streamline Lepton variable update process with ptr-vectors
srtee Jan 14, 2025
b27aa31
fix whitespace
srtee Jan 14, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
299 changes: 89 additions & 210 deletions src/LEPTON/fix_efield_lepton.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -179,37 +179,37 @@ void FixEfieldLepton::post_force(int vflag)
auto dphi_x = parsed.differentiate("x").createCompiledExpression();
auto dphi_y = parsed.differentiate("y").createCompiledExpression();
auto dphi_z = parsed.differentiate("z").createCompiledExpression();
std::vector<Lepton::CompiledExpression*> dphis = {&dphi_x, &dphi_y, &dphi_z};
std::array<Lepton::CompiledExpression*, 3> dphis = {&dphi_x, &dphi_y, &dphi_z};

// check if reference to x, y, z exist
const std::array<std::string, 3> variableNames = {"x", "y", "z"};
std::array<bool, 3> phi_has_ref = {true, true, true};
// array of vectors of ptrs to Lepton variable references
std::array<std::vector<double *>, 3> var_ref_ptrs{};

// fill ptr-vectors with Lepton refs as needed
const char* DIM_NAMES[] = {"x", "y", "z"};
if (atom->q_flag){
phi = parsed.createCompiledExpression();
for (size_t i = 0; i < 3; i++) {
for (size_t d = 0; d < 3; d++) {
try {
phi.getVariableReference(variableNames[i]);
}
catch (Lepton::Exception &) {
phi_has_ref[i] = false;
double *ptr = &(phi.getVariableReference(DIM_NAMES[d]));
var_ref_ptrs[d].push_back(ptr);
} catch (Lepton::Exception &) {
// do nothing
}
}
}
std::vector<std::array<bool, 3>> dphis_has_ref;

bool e_uniform = true;
for (auto &dphi : dphis) {
dphis_has_ref.push_back({false, false, false});
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < 3; j++)
for (size_t d = 0; d < 3; d++) {
try {
(*dphi).getVariableReference(variableNames[i]);
dphis_has_ref.back()[i] = true;
double *ptr = &((*dphis[j]).getVariableReference(DIM_NAMES[d]));
var_ref_ptrs[d].push_back(ptr);
e_uniform = false;
}
catch (Lepton::Exception &) {
// do nothing
}
}
}
if (!e_uniform && atom->mu_flag && h < 0) {
error->all(FLERR, "Fix {} requires keyword `step' for dipoles in a non-uniform electric field", style);
}
Expand All @@ -227,43 +227,45 @@ void FixEfieldLepton::post_force(int vflag)

double ex, ey, ez;
double fx, fy, fz;
double v[6], unwrap[3];
double v[6], unwrap[3], dstep[3];
double xf, yf, zf, xb, yb, zb;
double exf, eyf, ezf, exb, eyb, ezb;
double mu_norm, h_mu;

if (atom->q_flag && atom->mu_flag) {
double *q = atom->q;
double **mu = atom->mu;
double **t = atom->torque;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue;
domain->unmap(x[i], image[i], unwrap);

// evaluate e-field, used by q and mu
for (size_t j = 0; j < 3; j++) {
if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = unwrap[0];
if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = unwrap[1];
if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = unwrap[2];
double *q = atom->q;
double **mu = atom->mu;
double **t = atom->torque;

for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue;
fx = fy = fz = 0.0;
domain->unmap(x[i], image[i], unwrap);

// put unwrapped coords into Lepton variable refs
for (size_t d = 0; d < 3; d++) {
for (auto & var_ref_ptr : var_ref_ptrs[d]) {
*var_ref_ptr = unwrap[d];
}
ex = -dphi_x.evaluate();
ey = -dphi_y.evaluate();
ez = -dphi_z.evaluate();
}

if (phi_has_ref[0]) phi.getVariableReference("x") = unwrap[0];
if (phi_has_ref[1]) phi.getVariableReference("y") = unwrap[1];
if (phi_has_ref[2]) phi.getVariableReference("z") = unwrap[2];
// evaluate e-field, used by q and mu
ex = -dphi_x.evaluate();
ey = -dphi_y.evaluate();
ez = -dphi_z.evaluate();

// charges
// force = q E
// charges
// force = q E
if (atom->q_flag) {
fx = qe2f * q[i] * ex;
fy = qe2f * q[i] * ey;
fz = qe2f * q[i] * ez;
// potential energy = q phi
fsum[0] += qe2f * q[i] * phi.evaluate();
}

// dipoles
if (atom->mu_flag) {
// dipoles
mu_norm = sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]);
if (mu_norm > EPSILON) {
// torque = mu cross E
Expand All @@ -273,181 +275,58 @@ void FixEfieldLepton::post_force(int vflag)
// potential energy = - mu dot E
fsum[0] -= mue2e * (mu[i][0] * ex + mu[i][1] * ey + mu[i][2] * ez);

// force = (mu dot D) E
// force = (mu dot D) E for non-uniform E
// using central difference method
h_mu = h / mu_norm;

xf = unwrap[0] + h_mu * mu[i][0];
yf = unwrap[1] + h_mu * mu[i][1];
zf = unwrap[2] + h_mu * mu[i][2];
for (size_t j = 0; j < 3; j++) {
if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = xf;
if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = yf;
if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = zf;
}
exf = -dphi_x.evaluate();
eyf = -dphi_y.evaluate();
ezf = -dphi_z.evaluate();

xb = unwrap[0] - h_mu * mu[i][0];
yb = unwrap[1] - h_mu * mu[i][1];
zb = unwrap[2] - h_mu * mu[i][2];
for (size_t j = 0; j < 3; j++) {
if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = xb;
if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = yb;
if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = zb;
if (!e_uniform) {
h_mu = h / mu_norm;
dstep[0] = h_mu * mu[i][0];
dstep[1] = h_mu * mu[i][1];
dstep[2] = h_mu * mu[i][2];

// one step forwards, two steps back ;)
for (size_t d = 0; d < 3; d++) {
for (auto & var_ref_ptr : var_ref_ptrs[d]) {
*var_ref_ptr += dstep[d];
}
}

exf = -dphi_x.evaluate();
eyf = -dphi_y.evaluate();
ezf = -dphi_z.evaluate();

for (size_t d = 0; d < 3; d++) {
for (auto & var_ref_ptr : var_ref_ptrs[d]) {
*var_ref_ptr -= 2*dstep[d];
}
}

exb = -dphi_x.evaluate();
eyb = -dphi_y.evaluate();
ezb = -dphi_z.evaluate();

fx += qe2f * (exf - exb) / 2.0 / h_mu;
fy += qe2f * (eyf - eyb) / 2.0 / h_mu;
fz += qe2f * (ezf - ezb) / 2.0 / h_mu;
}
exb = -dphi_x.evaluate();
eyb = -dphi_y.evaluate();
ezb = -dphi_z.evaluate();

fx += qe2f * (exf - exb) / 2.0 / h_mu;
fy += qe2f * (eyf - eyb) / 2.0 / h_mu;
fz += qe2f * (ezf - ezb) / 2.0 / h_mu;
}

f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;

fsum[1] += fx;
fsum[2] += fy;
fsum[3] += fz;

if (evflag) {
v[0] = fx * unwrap[0];
v[1] = fy * unwrap[1];
v[2] = fz * unwrap[2];
v[3] = fx * unwrap[1];
v[4] = fx * unwrap[2];
v[5] = fy * unwrap[2];
v_tally(i, v);
}
}
}
} else if (atom->q_flag && !atom->mu_flag) {
double *q = atom->q;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue;
domain->unmap(x[i], image[i], unwrap);

for (size_t j = 0; j < 3; j++) {
if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = unwrap[0];
if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = unwrap[1];
if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = unwrap[2];
}
ex = -dphi_x.evaluate();
ey = -dphi_y.evaluate();
ez = -dphi_z.evaluate();

if (phi_has_ref[0]) phi.getVariableReference("x") = unwrap[0];
if (phi_has_ref[1]) phi.getVariableReference("y") = unwrap[1];
if (phi_has_ref[2]) phi.getVariableReference("z") = unwrap[2];

// force = q E
fx = qe2f * q[i] * ex;
fy = qe2f * q[i] * ey;
fz = qe2f * q[i] * ez;
// potential energy = q phi
fsum[0] += qe2f * q[i] * phi.evaluate();

f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;

fsum[1] += fx;
fsum[2] += fy;
fsum[3] += fz;

if (evflag) {
v[0] = fx * unwrap[0];
v[1] = fy * unwrap[1];
v[2] = fz * unwrap[2];
v[3] = fx * unwrap[1];
v[4] = fx * unwrap[2];
v[5] = fy * unwrap[2];
v_tally(i, v);
}
}
}
} else if (!atom->q_flag && atom->mu_flag) {
double **mu = atom->mu;
double **t = atom->torque;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue;

mu_norm = sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]);
if (mu_norm > EPSILON) continue;

domain->unmap(x[i], image[i], unwrap);

for (size_t j = 0; j < 3; j++) {
if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = unwrap[0];
if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = unwrap[1];
if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = unwrap[2];
}
ex = -dphi_x.evaluate();
ey = -dphi_y.evaluate();
ez = -dphi_z.evaluate();

// torque = mu cross E
t[i][0] += mue2e * (ez * mu[i][1] - ey * mu[i][2]);
t[i][1] += mue2e * (ex * mu[i][2] - ez * mu[i][0]);
t[i][2] += mue2e * (ey * mu[i][0] - ex * mu[i][1]);
// potential energy = - mu dot E
fsum[0] -= mue2e * (mu[i][0] * ex + mu[i][1] * ey + mu[i][2] * ez);

// force = (mu dot D) E
// using central difference method
h_mu = h / sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]);

xf = unwrap[0] + h_mu * mu[i][0];
yf = unwrap[1] + h_mu * mu[i][1];
zf = unwrap[2] + h_mu * mu[i][2];
for (size_t j = 0; j < 3; j++) {
if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = xf;
if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = yf;
if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = zf;
}
exf = -dphi_x.evaluate();
eyf = -dphi_y.evaluate();
ezf = -dphi_z.evaluate();

xb = unwrap[0] - h_mu * mu[i][0];
yb = unwrap[1] - h_mu * mu[i][1];
zb = unwrap[2] - h_mu * mu[i][2];
for (size_t j = 0; j < 3; j++) {
if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = xb;
if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = yb;
if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = zb;
}
exb = -dphi_x.evaluate();
eyb = -dphi_y.evaluate();
ezb = -dphi_z.evaluate();

fx = qe2f * (exf - exb) / 2.0 / h_mu;
fy = qe2f * (eyf - eyb) / 2.0 / h_mu;
fz = qe2f * (ezf - ezb) / 2.0 / h_mu;

f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;

fsum[1] += fx;
fsum[2] += fy;
fsum[3] += fz;

if (evflag) {
v[0] = fx * unwrap[0];
v[1] = fy * unwrap[1];
v[2] = fz * unwrap[2];
v[3] = fx * unwrap[1];
v[4] = fx * unwrap[2];
v[5] = fy * unwrap[2];
v_tally(i, v);
}
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;

fsum[1] += fx;
fsum[2] += fy;
fsum[3] += fz;

if (evflag) {
v[0] = fx * unwrap[0];
v[1] = fy * unwrap[1];
v[2] = fz * unwrap[2];
v[3] = fx * unwrap[1];
v[4] = fx * unwrap[2];
v[5] = fy * unwrap[2];
v_tally(i, v);
}
}
}
Expand Down