From 7b9f7be485d1a9f88d02f650a7fe1c3e39475f99 Mon Sep 17 00:00:00 2001 From: Shern Tee Date: Tue, 14 Jan 2025 14:47:33 +1000 Subject: [PATCH 1/4] replace vecs with arrays since size known at compile-time; make phi and dphi flagging uniform --- src/LEPTON/fix_efield_lepton.cpp | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/src/LEPTON/fix_efield_lepton.cpp b/src/LEPTON/fix_efield_lepton.cpp index 05912382a33..d9db3987635 100644 --- a/src/LEPTON/fix_efield_lepton.cpp +++ b/src/LEPTON/fix_efield_lepton.cpp @@ -179,37 +179,36 @@ 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 dphis = {&dphi_x, &dphi_y, &dphi_z}; + std::array dphis = {&dphi_x, &dphi_y, &dphi_z}; // check if reference to x, y, z exist - const std::array variableNames = {"x", "y", "z"}; - std::array phi_has_ref = {true, true, true}; + const char* DIM_NAMES[] = {"x", "y", "z"}; + std::array phi_has_ref{}; // zero-init if (atom->q_flag){ phi = parsed.createCompiledExpression(); for (size_t i = 0; i < 3; i++) { try { - phi.getVariableReference(variableNames[i]); + phi.getVariableReference(DIM_NAMES[i]); + phi_has_ref[i] = true; } catch (Lepton::Exception &) { - phi_has_ref[i] = false; + // do nothing } } } - std::vector> dphis_has_ref; + std::array, 3> dphis_has_ref{}; bool e_uniform = true; - for (auto &dphi : dphis) { - dphis_has_ref.push_back({false, false, false}); + for (size_t j = 0; j < 3; j++) for (size_t i = 0; i < 3; i++) { try { - (*dphi).getVariableReference(variableNames[i]); - dphis_has_ref.back()[i] = true; + (*dphis[j]).getVariableReference(DIM_NAMES[i]); + dphis_has_ref[j][i] = true; 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); } From 1f7533029be546591b05d79c3a2cfc7d310fdc66 Mon Sep 17 00:00:00 2001 From: Shern Tee Date: Tue, 14 Jan 2025 22:08:16 +1000 Subject: [PATCH 2/4] deduplicate force and torque calculations --- src/LEPTON/fix_efield_lepton.cpp | 275 +++++++++---------------------- 1 file changed, 78 insertions(+), 197 deletions(-) diff --git a/src/LEPTON/fix_efield_lepton.cpp b/src/LEPTON/fix_efield_lepton.cpp index d9db3987635..8ed27480709 100644 --- a/src/LEPTON/fix_efield_lepton.cpp +++ b/src/LEPTON/fix_efield_lepton.cpp @@ -231,38 +231,42 @@ void FixEfieldLepton::post_force(int vflag) 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]; - } - 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]; + 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); + + // 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]; + } + 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 + + 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]; 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 @@ -271,182 +275,59 @@ void FixEfieldLepton::post_force(int vflag) 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 + + // 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; + + 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; } - 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); } } } From 276b8d9c93967569e2085cee8edb57f957abf476 Mon Sep 17 00:00:00 2001 From: Shern Tee Date: Tue, 14 Jan 2025 22:31:42 +1000 Subject: [PATCH 3/4] streamline Lepton variable update process with ptr-vectors --- src/LEPTON/fix_efield_lepton.cpp | 75 ++++++++++++++++---------------- 1 file changed, 37 insertions(+), 38 deletions(-) diff --git a/src/LEPTON/fix_efield_lepton.cpp b/src/LEPTON/fix_efield_lepton.cpp index 8ed27480709..7bc2e433ed6 100644 --- a/src/LEPTON/fix_efield_lepton.cpp +++ b/src/LEPTON/fix_efield_lepton.cpp @@ -181,28 +181,29 @@ void FixEfieldLepton::post_force(int vflag) auto dphi_z = parsed.differentiate("z").createCompiledExpression(); std::array dphis = {&dphi_x, &dphi_y, &dphi_z}; - // check if reference to x, y, z exist + // array of vectors of ptrs to Lepton variable references + std::array, 3> var_ref_ptrs{}; + + // fill ptr-vectors with Lepton refs as needed const char* DIM_NAMES[] = {"x", "y", "z"}; - std::array phi_has_ref{}; // zero-init 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(DIM_NAMES[i]); - phi_has_ref[i] = true; - } - catch (Lepton::Exception &) { + double *ptr = &(phi.getVariableReference(DIM_NAMES[d])); + var_ref_ptrs[d].push_back(ptr); + } catch (Lepton::Exception &) { // do nothing } } } - std::array, 3> dphis_has_ref{}; + bool e_uniform = true; for (size_t j = 0; j < 3; j++) - for (size_t i = 0; i < 3; i++) { + for (size_t d = 0; d < 3; d++) { try { - (*dphis[j]).getVariableReference(DIM_NAMES[i]); - dphis_has_ref[j][i] = true; + double *ptr = &((*dphis[j]).getVariableReference(DIM_NAMES[d])); + var_ref_ptrs[d].push_back(ptr); e_uniform = false; } catch (Lepton::Exception &) { @@ -226,7 +227,7 @@ 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; @@ -241,12 +242,14 @@ void FixEfieldLepton::post_force(int vflag) fx = fy = fz = 0.0; 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]; + // 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]; + } } + + // evaluate e-field, used by q and mu ex = -dphi_x.evaluate(); ey = -dphi_y.evaluate(); ez = -dphi_z.evaluate(); @@ -258,10 +261,6 @@ void FixEfieldLepton::post_force(int vflag) fy = qe2f * q[i] * ey; fz = qe2f * q[i] * ez; // potential energy = q phi - - 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]; fsum[0] += qe2f * q[i] * phi.evaluate(); } @@ -280,27 +279,27 @@ void FixEfieldLepton::post_force(int vflag) // using central difference method if (!e_uniform) { 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; + 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(); + + 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; + + 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(); From b27aa31baa859aa47e9499144b2eaeb0846385d4 Mon Sep 17 00:00:00 2001 From: Shern Tee Date: Tue, 14 Jan 2025 22:35:56 +1000 Subject: [PATCH 4/4] fix whitespace --- src/LEPTON/fix_efield_lepton.cpp | 34 ++++++++++++++++---------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/LEPTON/fix_efield_lepton.cpp b/src/LEPTON/fix_efield_lepton.cpp index 7bc2e433ed6..a055c2b15ae 100644 --- a/src/LEPTON/fix_efield_lepton.cpp +++ b/src/LEPTON/fix_efield_lepton.cpp @@ -183,7 +183,7 @@ void FixEfieldLepton::post_force(int vflag) // array of vectors of ptrs to Lepton variable references std::array, 3> var_ref_ptrs{}; - + // fill ptr-vectors with Lepton refs as needed const char* DIM_NAMES[] = {"x", "y", "z"}; if (atom->q_flag){ @@ -197,7 +197,7 @@ void FixEfieldLepton::post_force(int vflag) } } } - + bool e_uniform = true; for (size_t j = 0; j < 3; j++) for (size_t d = 0; d < 3; d++) { @@ -274,41 +274,41 @@ void FixEfieldLepton::post_force(int vflag) 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 for non-uniform E // using central difference method - if (!e_uniform) { + 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++) { + 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(); + + exf = -dphi_x.evaluate(); eyf = -dphi_y.evaluate(); ezf = -dphi_z.evaluate(); - - for (size_t d = 0; d < 3; d++) { + + 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; } - } + } } f[i][0] += fx;