espressomd / espresso

The ESPResSo package
https://espressomd.org
GNU General Public License v3.0
224 stars 183 forks source link

Separate generic properties from interaction potentials #3060

Open fweik opened 5 years ago

fweik commented 5 years ago

The actual interaction potentials should be separated from generic properties like cutoff, offset and shift. This is independent from each other, would give all potentials shifts and cutoff and would reduce code duplication. We also had bugs in the past because potentials did not correctly calculate the cutoff if there was a shift. Here is a rough draft of what I have in mind:

template<class Potential>
struct CetralPotential : Potential {
   double shift, offset, cutoff;

   double energy(double r) {
       auto const r_eff = r + offset;
      if(r_eff <= cutoff) 
        return Potential::U(r) + shift;
      else
        return 0.;
   }
};

struct Harmonic {
  double k;
   double U(double r) {
       return -k*r*r;
   }
};

struct IA_Parameters {
  [...]
  CentralPotential<Harmonic> harmonic;
  [...]
}

This also allows for easier reuse of the potential functions (observe how the potential knows nothing about how it is used). If all the potentials have the same interface, they can also be used in a homogeneous way, e.g. we can have

using IA_Parameters = std::tuple<Potentials...>;
IA_Parameters ia_params;

double energy = 0.;
Utils::for_each(ia_params, [r, energy](auto const& ia) {
  energy += ia.energy(r);
});  

Some more thought should go into this I think, but I think something along these lines is a good starting point.

KaiSzuttor commented 5 years ago
auto const energy = std::accumulate(ia_params.begin(), ia_params.end(), 0.0, [r](double e, CentralPotential const &p) { return e + p.energy(r); }

would also work?

fweik commented 5 years ago

No. (it says Utils::for_each not std::for_each for a reason). Also discussions of implementation details make more sense if there is an actual implementation.

KaiSzuttor commented 5 years ago

it was just for curiosity not for productivity

fweik commented 5 years ago

STL algorithms only work on homogeneous containers, in your code CentralPotential is a template, not a type, it would already fail there.

KaiSzuttor commented 5 years ago
auto const energy = boost::fusion::accumulate(ia_params, 0.0, [r](double e, auto const &p) { return e + p.energy(r); }

should work

KaiSzuttor commented 5 years ago

side note: it should be auto const r_eff = r - offset

KaiSzuttor commented 5 years ago

what about using something like this for the interactions container?

struct CantorPairing {
    std::size_t cantor_pairing(int k1, int k2) {
      return static_cast<std::size_t>(0.5 * (k1 + k2) * (k1 + k2 + 1) + k2);
    }

    std::size_t operator()(std::pair<int, int> p) {
      // make it symmetric
      if (p.second > p.first) {
        return cantor_pairing(p.first, p.second);
      }
      return cantor_pairing(p.second, p.first);
    }
};

std::unordered_map<std::pair<int, int>, IA_parameters, CantorPairing> interactions;
fweik commented 5 years ago

What are you trying to do here? Not sure we are on the same page...

KaiSzuttor commented 5 years ago

just experimenting with the mapping of type pairs to interactions

fweik commented 5 years ago

This is not what this is about.

On Thu, Aug 8, 2019, 18:04 Kai Szuttor notifications@github.com wrote:

just experimenting with the mapping of type pairs to interactions

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/espressomd/espresso/issues/3060?email_source=notifications&email_token=AAG2FX5VAXFSUYS7J6GVC7TQDQ7YVA5CNFSM4IJ7PFA2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD34DFUQ#issuecomment-519582418, or mute the thread https://github.com/notifications/unsubscribe-auth/AAG2FXZHQABMU6O42Y3ALJLQDQ7YVANCNFSM4IJ7PFAQ .

KaiSzuttor commented 5 years ago

i see

Am 08.08.2019 um 18:06 schrieb Florian Weik notifications@github.com:

This is not what this is about.

On Thu, Aug 8, 2019, 18:04 Kai Szuttor notifications@github.com wrote:

just experimenting with the mapping of type pairs to interactions

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/espressomd/espresso/issues/3060?email_source=notifications&email_token=AAG2FX5VAXFSUYS7J6GVC7TQDQ7YVA5CNFSM4IJ7PFA2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD34DFUQ#issuecomment-519582418, or mute the thread https://github.com/notifications/unsubscribe-auth/AAG2FXZHQABMU6O42Y3ALJLQDQ7YVANCNFSM4IJ7PFAQ .

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or mute the thread.

KaiSzuttor commented 5 years ago

could you please explain what data structure you had in mind for the potentials? you used a struct and tuple for IA_parameters in your snippets got it

fweik commented 5 years ago

Just to spell it out: this is about multiple different interaction e.g. per pair or bond. This is independent of how those interactions are assigned to the interaction pairs, e.g. by type for which there would be an separate data structure, containing multiple such tuples as discussed here. Sorry if that wasn't clear.