Skip to content

Commit

Permalink
Refactor multipole distribution
Browse files Browse the repository at this point in the history
  • Loading branch information
@mlund
mlund committed Nov 30, 2023
1 parent 06f1e83 commit 00ce667
Show file tree
Hide file tree
Showing 2 changed files with 50 additions and 47 deletions.
88 changes: 46 additions & 42 deletions src/analysis.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#include "analysis.h"
#include "auxiliary.h"
#include "core.h"
#include "move.h"
#include "energy.h"
Expand Down Expand Up @@ -1472,61 +1473,66 @@ void XTCtraj::_sample() {

// =============== MultipoleDistribution ===============

double MultipoleDistribution::g2g(const Space::GroupType& group1, const Space::GroupType& group2) {
double MultipoleDistribution::groupGroupExactEnergy(const Space::GroupType& group1,
const Space::GroupType& group2) const {
double energy = 0.0;
for (const auto& particle_i : group1) {
for (const auto& particle_j : group2) {
energy += particle_i.charge * particle_j.charge / spc.geometry.vdist(particle_i.pos, particle_j.pos).norm();
for (const auto& a : group1) {
for (const auto& b : group2) {
energy += a.charge * b.charge / spc.geometry.vdist(a.pos, b.pos).norm();
}
}
return energy;
}

/**
* @note `fmt` is currently included w. spdlog but has been accepted into c++20.
*/
void MultipoleDistribution::save() const {
if (number_of_samples > 0) {
if (std::ofstream stream(MPI::prefix + filename.c_str()); stream) {
stream << "# Multipolar energies (kT/lB)\n"
<< fmt::format("# {:>8}{:>10}{:>10}{:>10}{:>10}{:>10}{:>10}{:>10}\n", "R", "exact", "tot", "ii",
"id", "dd", "iq", "mucorr");
for (auto [r_bin, energy] : mean_energy) {
const auto distance = r_bin * dr;
const auto u_tot = energy.ion_ion.avg() + energy.ion_dipole.avg() + energy.dipole_dipole.avg() +
energy.ion_quadrupole.avg();
stream << fmt::format("{:10.4f}{:10.4f}{:10.4f}{:10.4f}{:10.4f}{:10.4f}{:10.4f}{:10.4f}\n", distance,
energy.exact.avg(), u_tot, energy.ion_ion.avg(), energy.ion_dipole.avg(),
energy.dipole_dipole.avg(), energy.ion_quadrupole.avg(),
energy.dipole_dipole_correlation.avg());
}
}
}
/// Calculate multipole energies between two groups and update averages at given separation
void MultipoleDistribution::sampleGroupGroup(const Space::GroupType& group1, const Space::GroupType& group2) {
const auto a = Faunus::toMultipole(group1, spc.geometry.getBoundaryFunc());
const auto b = Faunus::toMultipole(group2, spc.geometry.getBoundaryFunc());
const auto distance = spc.geometry.vdist(group1.mass_center, group2.mass_center);
auto& data = mean_energy[to_bin(distance.norm(), dr)];
data.exact += groupGroupExactEnergy(group1, group2);
data.ion_ion += a.charge * b.charge / distance.norm();
data.ion_dipole +=
q2mu(a.charge * b.getExt().mulen, b.getExt().mu, b.charge * a.getExt().mulen, a.getExt().mu, distance);
data.dipole_dipole += mu2mu(a.getExt().mu, b.getExt().mu, a.getExt().mulen * b.getExt().mulen, distance);
data.ion_quadrupole += q2quad(a.charge, b.getExt().Q, b.charge, a.getExt().Q, distance);
data.dipole_dipole_correlation += a.getExt().mu.dot(b.getExt().mu);
}

void MultipoleDistribution::_sample() {
for (auto& group1 : spc.findMolecules(ids[0])) {
// find active molecules
for (auto& group2 : spc.findMolecules(ids[1])) {
// find active molecules
// loop over active molecules
for (const auto& group1 : spc.findMolecules(ids.at(0))) {
for (const auto& group2 : spc.findMolecules(ids.at(1))) {
if (&group1 != &group2) {
const auto a = Faunus::toMultipole(group1, spc.geometry.getBoundaryFunc());
const auto b = Faunus::toMultipole(group2, spc.geometry.getBoundaryFunc());
const auto distance = spc.geometry.vdist(group1.mass_center, group2.mass_center);
auto& data = mean_energy[to_bin(distance.norm(), dr)];
data.exact += g2g(group1, group2);
data.ion_ion += a.charge * b.charge / distance.norm();
data.ion_dipole += q2mu(a.charge * b.getExt().mulen, b.getExt().mu, b.charge * a.getExt().mulen,
a.getExt().mu, distance);
data.dipole_dipole +=
mu2mu(a.getExt().mu, b.getExt().mu, a.getExt().mulen * b.getExt().mulen, distance);
data.ion_quadrupole += q2quad(a.charge, b.getExt().Q, b.charge, a.getExt().Q, distance);
data.dipole_dipole_correlation += a.getExt().mu.dot(b.getExt().mu);
sampleGroupGroup(group1, group2);
}
}
}
}

/**
* @note `fmt` is currently included w. spdlog but has been accepted into c++20.
*/
void MultipoleDistribution::_to_disk() {
if (number_of_samples == 0) {
return;
}
if (std::ofstream stream(MPI::prefix + filename.c_str()); stream) {
stream << "# Multipolar energies (kT/lB)\n"
<< fmt::format("# {:>8}{:>10}{:>10}{:>10}{:>10}{:>10}{:>10}{:>10}\n", "R", "exact", "tot", "ii", "id",
"dd", "iq", "mucorr");
for (auto [r_bin, energy] : mean_energy) {
const auto distance = r_bin * dr;
const auto u_tot = energy.ion_ion.avg() + energy.ion_dipole.avg() + energy.dipole_dipole.avg() +
energy.ion_quadrupole.avg();
stream << fmt::format("{:10.4f}{:10.4f}{:10.4f}{:10.4f}{:10.4f}{:10.4f}{:10.4f}{:10.4f}\n", distance,
energy.exact.avg(), u_tot, energy.ion_ion.avg(), energy.ion_dipole.avg(),
energy.dipole_dipole.avg(), energy.ion_quadrupole.avg(),
energy.dipole_dipole_correlation.avg());
}
}
}

void MultipoleDistribution::_to_json(json& j) const { j = {{"molecules", names}, {"file", filename}, {"dr", dr}}; }

MultipoleDistribution::MultipoleDistribution(const json& j, const Space& spc)
Expand All @@ -1541,8 +1547,6 @@ MultipoleDistribution::MultipoleDistribution(const json& j, const Space& spc)
}
}

void MultipoleDistribution::_to_disk() { save(); }

// =============== AtomInertia ===============

void AtomInertia::_to_json(json& j) const {
Expand Down
9 changes: 4 additions & 5 deletions src/analysis.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -733,7 +733,7 @@ class VirtualTranslate : public PerturbationAnalysis {
class MultipoleDistribution : public Analysis {
struct Data {
using average_type = Average<double>;
average_type exact;
average_type exact; //!< Exact electrostatic energy
average_type ion_ion;
average_type ion_dipole;
average_type ion_quadrupole;
Expand All @@ -746,16 +746,15 @@ class MultipoleDistribution : public Analysis {
std::string filename; //!< output file name
double dr = 0.0; //!< distance resolution

double g2g(const Space::GroupType& group1,
const Space::GroupType& group2); //<! exact ion-ion energy between particles
void save() const; //!< save to disk
void _sample() override;
void _to_json(json& j) const override;
void _to_disk() override;
void sampleGroupGroup(const Space::GroupType& group1, const Space::GroupType& group2);
double groupGroupExactEnergy(const Space::GroupType& group1,
const Space::GroupType& group2) const; //<! exact ion-ion energy between particles

public:
MultipoleDistribution(const json& j, const Space& spc);

}; // end of multipole distribution

/**
Expand Down

0 comments on commit 00ce667

Please sign in to comment.