Skip to content

Commit

Permalink
core: elc: removed custom inverse box dimensions
Browse files Browse the repository at this point in the history
  • Loading branch information
@reinaual
reinaual committed May 18, 2021
1 parent 354f540 commit 05f5a06
Showing 1 changed file with 73 additions and 56 deletions.
129 changes: 73 additions & 56 deletions src/core/electrostatics_magnetostatics/elc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,11 +50,6 @@
#include <cstddef>
#include <vector>

/** \name Inverse box dimensions */
/**@{*/
static double ux, uy, uz;
/**@}*/

ELC_struct elc_params = {1e100, 10, 1, 0, true, true, false, 1,
1, false, 0, 0, 0, 0, 0.0};

Expand Down Expand Up @@ -107,12 +102,6 @@ static double dipole_energy(const ParticleRange &particles);
static double z_energy(const ParticleRange &particles);
static void add_z_force(const ParticleRange &particles);

void ELC_setup_constants() {
ux = 1 / box_geo.length()[0];
uy = 1 / box_geo.length()[1];
uz = 1 / box_geo.length()[2];
}

/**
* @brief Calculate cached sin/cos values for one direction.
*
Expand Down Expand Up @@ -215,7 +204,9 @@ inline void check_gap_elc(const Particle &p) {
* See @cite yeh99a.
*/
static void add_dipole_force(const ParticleRange &particles) {
double const pref = coulomb.prefactor * 4 * Utils::pi() * ux * uy * uz;
double const pref = coulomb.prefactor * 4 * Utils::pi() *
box_geo.length_inv()[0] * box_geo.length_inv()[1] *
box_geo.length_inv()[2];
constexpr std::size_t size = 3;

auto local_particles = particles;
Expand Down Expand Up @@ -251,7 +242,7 @@ static void add_dipole_force(const ParticleRange &particles) {
}

gblcblk[0] *= pref;
gblcblk[1] *= pref / elc_params.h / uz;
gblcblk[1] *= pref / elc_params.h * box_geo.length()[2];
gblcblk[2] *= pref;

distribute(size);
Expand Down Expand Up @@ -280,7 +271,9 @@ static void add_dipole_force(const ParticleRange &particles) {
* See @cite yeh99a.
*/
static double dipole_energy(const ParticleRange &particles) {
double const pref = coulomb.prefactor * 2 * Utils::pi() * ux * uy * uz;
double const pref = coulomb.prefactor * 2 * Utils::pi() *
box_geo.length_inv()[0] * box_geo.length_inv()[1] *
box_geo.length_inv()[2];
constexpr std::size_t size = 7;
/* for nonneutral systems, this shift gives the background contribution
(rsp. for this shift, the DM of the background is zero) */
Expand Down Expand Up @@ -336,7 +329,8 @@ static double dipole_energy(const ParticleRange &particles) {
if (elc_params.dielectric_contrast_on) {
if (elc_params.const_pot) {
// zero potential difference contribution
energy += pref / elc_params.h / uz * Utils::sqr(gblcblk[6]);
energy +=
pref / elc_params.h * box_geo.length()[2] * Utils::sqr(gblcblk[6]);
// external potential shift contribution
energy -= 2 * elc_params.pot_diff / elc_params.h * gblcblk[6];
}
Expand Down Expand Up @@ -374,7 +368,8 @@ inline double image_sum_t(double q, double z) {

/*****************************************************************/
static double z_energy(const ParticleRange &particles) {
double const pref = coulomb.prefactor * 2 * Utils::pi() * ux * uy;
double const pref = coulomb.prefactor * 2 * Utils::pi() *
box_geo.length_inv()[0] * box_geo.length_inv()[1];
constexpr std::size_t size = 4;

/* for nonneutral systems, this shift gives the background contribution
Expand Down Expand Up @@ -452,7 +447,8 @@ static double z_energy(const ParticleRange &particles) {

/*****************************************************************/
static void add_z_force(const ParticleRange &particles) {
double const pref = coulomb.prefactor * 2 * Utils::pi() * ux * uy;
double const pref = coulomb.prefactor * 2 * Utils::pi() *
box_geo.length_inv()[0] * box_geo.length_inv()[1];
constexpr std::size_t size = 1;

if (elc_params.dielectric_contrast_on) {
Expand Down Expand Up @@ -512,7 +508,8 @@ template <PoQ axis>
void setup_PoQ(std::size_t index, double omega,
const ParticleRange &particles) {
assert(index >= 1);
double const pref_di = coulomb.prefactor * 4 * Utils::pi() * ux * uy;
double const pref_di = coulomb.prefactor * 4 * Utils::pi() *
box_geo.length_inv()[0] * box_geo.length_inv()[1];
double const pref = -pref_di / expm1(omega * box_geo.length()[2]);
constexpr std::size_t size = 4;
double lclimgebot[4], lclimgetop[4], lclimge[4];
Expand Down Expand Up @@ -659,7 +656,8 @@ static void setup_PQ(std::size_t index_p, std::size_t index_q, double omega,
const ParticleRange &particles) {
assert(index_p >= 1);
assert(index_q >= 1);
double const pref_di = coulomb.prefactor * 8 * Utils::pi() * ux * uy;
double const pref_di = coulomb.prefactor * 8 * Utils::pi() *
box_geo.length_inv()[0] * box_geo.length_inv()[1];
double const pref = -pref_di / expm1(omega * box_geo.length()[2]);
constexpr std::size_t size = 8;
double lclimgebot[8], lclimgetop[8], lclimge[8];
Expand Down Expand Up @@ -789,8 +787,10 @@ static void setup_PQ(std::size_t index_p, std::size_t index_q, double omega,
static void add_PQ_force(std::size_t index_p, std::size_t index_q, double omega,
const ParticleRange &particles) {
constexpr double c_2pi = 2 * Utils::pi();
double const pref_x = c_2pi * ux * static_cast<double>(index_p) / omega;
double const pref_y = c_2pi * uy * static_cast<double>(index_q) / omega;
double const pref_x =
c_2pi * box_geo.length_inv()[0] * static_cast<double>(index_p) / omega;
double const pref_y =
c_2pi * box_geo.length_inv()[1] * static_cast<double>(index_q) / omega;
constexpr std::size_t size = 8;

std::size_t ic = 0;
Expand Down Expand Up @@ -847,45 +847,54 @@ static double PQ_energy(double omega, std::size_t n_part) {

void ELC_add_force(const ParticleRange &particles) {
constexpr double c_2pi = 2 * Utils::pi();
auto const n_scxcache = std::size_t(ceil(elc_params.far_cut / ux) + 1);
auto const n_scycache = std::size_t(ceil(elc_params.far_cut / uy) + 1);
auto const n_scxcache =
std::size_t(ceil(elc_params.far_cut * box_geo.length()[0]) + 1);
auto const n_scycache =
std::size_t(ceil(elc_params.far_cut * box_geo.length()[1]) + 1);

prepare_sc_cache(particles, n_scxcache, ux, n_scycache, uy);
prepare_sc_cache(particles, n_scxcache, box_geo.length_inv()[0], n_scycache,
box_geo.length_inv()[1]);
partblk.resize(particles.size() * 8);

add_dipole_force(particles);
add_z_force(particles);

/* the second condition is just for the case of numerical accident */
for (std::size_t p = 1;
ux * static_cast<double>(p - 1) < elc_params.far_cut && p <= n_scxcache;
for (std::size_t p = 1; box_geo.length_inv()[0] * static_cast<double>(p - 1) <
elc_params.far_cut &&
p <= n_scxcache;
p++) {
auto const omega = c_2pi * ux * static_cast<double>(p);
auto const omega = c_2pi * box_geo.length_inv()[0] * static_cast<double>(p);
setup_PoQ<PoQ::P>(p, omega, particles);
distribute(4);
add_PoQ_force<PoQ::P>(particles);
}

for (std::size_t q = 1;
uy * static_cast<double>(q - 1) < elc_params.far_cut && q <= n_scycache;
for (std::size_t q = 1; box_geo.length_inv()[1] * static_cast<double>(q - 1) <
elc_params.far_cut &&
q <= n_scycache;
q++) {
auto const omega = c_2pi * uy * static_cast<double>(q);
auto const omega = c_2pi * box_geo.length_inv()[1] * static_cast<double>(q);
setup_PoQ<PoQ::Q>(q, omega, particles);
distribute(4);
add_PoQ_force<PoQ::Q>(particles);
}

for (std::size_t p = 1;
ux * static_cast<double>(p - 1) < elc_params.far_cut && p <= n_scxcache;
for (std::size_t p = 1; box_geo.length_inv()[0] * static_cast<double>(p - 1) <
elc_params.far_cut &&
p <= n_scxcache;
p++) {
for (std::size_t q = 1;
Utils::sqr(ux * static_cast<double>(p - 1)) +
Utils::sqr(uy * static_cast<double>(q - 1)) <
Utils::sqr(box_geo.length_inv()[0] * static_cast<double>(p - 1)) +
Utils::sqr(box_geo.length_inv()[1] *
static_cast<double>(q - 1)) <
elc_params.far_cut2 &&
q <= n_scycache;
q++) {
auto const omega = c_2pi * sqrt(Utils::sqr(ux * static_cast<double>(p)) +
Utils::sqr(uy * static_cast<double>(q)));
auto const omega =
c_2pi *
sqrt(Utils::sqr(box_geo.length_inv()[0] * static_cast<double>(p)) +
Utils::sqr(box_geo.length_inv()[1] * static_cast<double>(q)));
setup_PQ(p, q, omega, particles);
distribute(8);
add_PQ_force(p, q, omega, particles);
Expand All @@ -898,43 +907,52 @@ double ELC_energy(const ParticleRange &particles) {
auto energy = dipole_energy(particles);
energy += z_energy(particles);

auto const n_scxcache = std::size_t(ceil(elc_params.far_cut / ux) + 1);
auto const n_scycache = std::size_t(ceil(elc_params.far_cut / uy) + 1);
prepare_sc_cache(particles, n_scxcache, ux, n_scycache, uy);
auto const n_scxcache =
std::size_t(ceil(elc_params.far_cut * box_geo.length()[0]) + 1);
auto const n_scycache =
std::size_t(ceil(elc_params.far_cut * box_geo.length()[1]) + 1);
prepare_sc_cache(particles, n_scxcache, box_geo.length_inv()[0], n_scycache,
box_geo.length_inv()[1]);

auto const n_localpart = particles.size();
partblk.resize(n_localpart * 8);

/* the second condition is just for the case of numerical accident */
for (std::size_t p = 1;
ux * static_cast<double>(p - 1) < elc_params.far_cut && p <= n_scxcache;
for (std::size_t p = 1; box_geo.length_inv()[0] * static_cast<double>(p - 1) <
elc_params.far_cut &&
p <= n_scxcache;
p++) {
auto const omega = c_2pi * ux * static_cast<double>(p);
auto const omega = c_2pi * box_geo.length_inv()[0] * static_cast<double>(p);
setup_PoQ<PoQ::P>(p, omega, particles);
distribute(4);
energy += PoQ_energy(omega, n_localpart);
}

for (std::size_t q = 1;
uy * static_cast<double>(q - 1) < elc_params.far_cut && q <= n_scycache;
for (std::size_t q = 1; box_geo.length_inv()[1] * static_cast<double>(q - 1) <
elc_params.far_cut &&
q <= n_scycache;
q++) {
auto const omega = c_2pi * uy * static_cast<double>(q);
auto const omega = c_2pi * box_geo.length_inv()[1] * static_cast<double>(q);
setup_PoQ<PoQ::Q>(q, omega, particles);
distribute(4);
energy += PoQ_energy(omega, n_localpart);
}

for (std::size_t p = 1;
ux * static_cast<double>(p - 1) < elc_params.far_cut && p <= n_scxcache;
for (std::size_t p = 1; box_geo.length_inv()[0] * static_cast<double>(p - 1) <
elc_params.far_cut &&
p <= n_scxcache;
p++) {
for (std::size_t q = 1;
Utils::sqr(ux * static_cast<double>(p - 1)) +
Utils::sqr(uy * static_cast<double>(q - 1)) <
Utils::sqr(box_geo.length_inv()[0] * static_cast<double>(p - 1)) +
Utils::sqr(box_geo.length_inv()[1] *
static_cast<double>(q - 1)) <
elc_params.far_cut2 &&
q <= n_scycache;
q++) {
auto const omega = c_2pi * sqrt(Utils::sqr(ux * static_cast<double>(p)) +
Utils::sqr(uy * static_cast<double>(q)));
auto const omega =
c_2pi *
sqrt(Utils::sqr(box_geo.length_inv()[0] * static_cast<double>(p)) +
Utils::sqr(box_geo.length_inv()[1] * static_cast<double>(q)));
setup_PQ(p, q, omega, particles);
distribute(8);
energy += PQ_energy(omega, n_localpart);
Expand All @@ -949,7 +967,8 @@ double ELC_tune_far_cut(ELC_struct const &params) {
constexpr auto maximal_far_cut = 50.;
double const h = params.h;
double lz = box_geo.length()[2];
double const min_inv_boxl = std::min(ux, uy);
double const min_inv_boxl =
std::min(box_geo.length_inv()[0], box_geo.length_inv()[1]);

if (params.dielectric_contrast_on) {
// adjust lz according to dielectric layer method
Expand All @@ -965,7 +984,8 @@ double ELC_tune_far_cut(ELC_struct const &params) {
do {
const auto prefactor = 2 * Utils::pi() * far_cut;

const auto sum = prefactor + 2 * (ux + uy);
const auto sum =
prefactor + 2 * (box_geo.length_inv()[0] + box_geo.length_inv()[1]);
const auto den = -expm1(-prefactor * lz);
const auto num1 = exp(prefactor * (h - lz));
const auto num2 = exp(-prefactor * (h + lz));
Expand Down Expand Up @@ -1018,8 +1038,6 @@ void ELC_sanity_checks(ELC_struct const &params) {
}

void ELC_init() {

ELC_setup_constants();
elc_params.h = box_geo.length()[2] - elc_params.gap_size;

if (elc_params.dielectric_contrast_on) {
Expand Down Expand Up @@ -1107,7 +1125,6 @@ void ELC_set_params(double maxPWerror, double gap_size, double far_cut,

ELC_sanity_checks(new_elc_params);

ELC_setup_constants();
if (new_elc_params.far_calculated) {
new_elc_params.far_cut = ELC_tune_far_cut(new_elc_params);
}
Expand Down

0 comments on commit 05f5a06

Please sign in to comment.