Skip to content

Instantly share code, notes, and snippets.

@mlund

mlund/tabulator.cpp

Last active December 14, 2015 22:59
Show Gist options
  • Save mlund/5162296 to your computer and use it in GitHub Desktop.
Save mlund/5162296 to your computer and use it in GitHub Desktop.
Download ZIP
Possible design for tabulated pair potential. Require C++11 flags, for example "g++ -std=c++11"
Raw
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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment