Possible design for tabulated pair potential. Require C++11 flags, for example "g++ -std=c++11"

tabulator.cpp

#include <iostream>
#include <functional>
#include <map>
#include <cmath>
#include <memory>

/* ordered pair */
template<class T>
struct opair : public std::pair<T,T> {
  typedef std::pair<T,T> base;
  opair() {}
  opair(const T &a, const T &b) : base(a,b) {
    if (a>b)
      std::swap(base::first,base::second);
  }
  bool find(const T &i) const {
    assert(base::first<=base::second);
    if (i!=base::first)
      if (i!=base::second)
        return false;
    return true;
  }
};

struct particle {
  double x,y,z,q,r; // etc...
  int id;
};

struct PairPotentialBase {
  virtual double operator()(const particle&, const particle&, double)=0;
};

struct coulomb : PairPotentialBase {
  double operator()(const particle &a, const particle &b, double r2) {
    return a.q*b.q/std::sqrt(r2);
  }
};

struct lennardjones : PairPotentialBase {
  double operator()(const particle &a, const particle &b, double r2) {
    double x=r2*r2*r2;
    return 1/(x*x) - 1/x;
  }
};

/* base class for all tabulators - no dependencies */
template<typename T=double>
class tabulatorbase {
  protected:
    T utol, ftol;
  public:
    struct data {
      T rmin2, rmax2; // useful to save these with table
      T a,b,c;
    };
    void setTolerance(T utol, T ftol=-1);
    void setRange(T rmin, T rmax) {}
};

/* tabulator 1 */
template<typename T=double>
struct tabulator : public tabulatorbase<T> {
  typedef tabulatorbase<T> base; // for convenience, only

  T eval(const typename base::data&, T r2) const {return 0;}  // gets tabulated value at r2

  typename base::data generate(std::function<T(T)> f) {
    std::cout << f(2) << "\n";
  }
};

/* THIS ALREADY EXISTS IN FAUNUS */
template<typename Tdefault, typename Tpair=opair<int>,
  typename Tmap=std::map<Tpair, std::shared_ptr<PairPotentialBase> > >
  class PotentialMap : public Tdefault {
    protected:
      Tmap m;
    public:
      template<class Tpairpot>
        void add(int id1, int id2, Tpairpot pot) {
          m[Tpair(id1,id2)] = std::shared_ptr<Tpairpot>(new Tpairpot(pot));
        }

      template<class Tparticle, class Tdist>
        double operator()(const Tparticle &a, const Tparticle &b, Tdist r2) {
          auto i=m.find( Tpair(a.id,b.id) );
          if (i!=m.end())
            return i->second->operator()(a,b,r2);
          return Tdefault::operator()(a,b,r2);
        }
  };

template<typename Tdefault, typename Ttabulator=tabulator<double>, typename Tpair=opair<int> >
class PotentialMapTabulated : public PotentialMap<Tdefault,Tpair> {
  private:
    typedef PotentialMap<Tdefault,Tpair> base;
    Ttabulator tab;
    std::map<Tpair, typename Ttabulator::data> mtab;
  public:
    PotentialMapTabulated() {
      // read from input file instead
      tab.setRange(2,100);
      particle a,b;
    }
    template<class Tparticle>
      double operator()(const Tparticle &a, const Tparticle &b, double r2) {
        auto ab = Tpair(a.id,b.id);
        auto it=mtab.find(ab);
        if (it!=mtab.end()) {
          if (r2<it->second.rmax2)
            if (r2>it->second.rmin2)
              return tab.eval(it->second, r2);
          return base::m[ab]->operator()(a,b,r2); // fall back to original
        }
        return Tdefault::operator()(a,b,r2); // fall back to default
      }
    template<class Tparticle, class Tpairpot>
      void add(Tparticle a, Tparticle b, Tpairpot pot) {
        base::add(a.id,b.id,pot);
        std::function<double(double)> f = [=](double r2) {return Tpairpot(pot)(a,b,r2);};
        mtab[ Tpair(a.id,b.id) ] = tab.generate(f);
      }
};

template<typename Tpairpot, typename Ttabulator=tabulator<double> >
class Tabulate : public Tpairpot {
  private:
    Ttabulator tab;
    typedef opair<int> Tpair;
    std::map<Tpair, typename Ttabulator::data > m;
  public:
    double operator()(const particle &a, const particle &b, double r2) {
      Tpair ab(a.id,b.id);
      auto it=m.find(ab);
      if (it!=m.end())
        return tab.eval(it->second, r2);
      m[ab] = tab.generate( [=](double x){ return Tpairpot(*this)(a,b,x); } );
      return (*this)(a,b,r2);
    }
};

int main() {

  particle a,b;
  tabulator<double> t;

  t.setRange(3.5, 200);
  t.generate( [=](double r2){ return coulomb()(a,b,r2); } );
  t.generate( [=](double r2){ return lennardjones()(a,b,r2); } );
  t.generate( acosh );
  t.generate( [](double x){ return 1/x; } );

  PotentialMapTabulated<coulomb> tt;
  tt.add(a,b,lennardjones());
  tt.add(b,b,coulomb());
  tt(a,b,20);

  //Tabulate<coulomb> coulombtab;

}