pele
Python energy landscape explorer
/home/js850/projects/pele/source/pele/distance.h
00001 #ifndef _PELE_DISTANCE_H
00002 #define _PELE_DISTANCE_H
00003 
00004 #include <cmath>
00005 #include <stdexcept>
00006 #include <type_traits>
00007 
00008 #include "array.h"
00009 
00016 // round is missing in visual studio
00017 #ifdef _MSC_VER
00018     inline double round(double r) {
00019         return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
00020     }
00021 #endif
00022 
00038 namespace pele {
00039 
00043 template<size_t IDX>
00044 struct meta_dist {
00045     static void f(double * const r_ij, double const * const r1,
00046             double const * const r2)
00047     {
00048         const static size_t k = IDX - 1;
00049         r_ij[k] = r1[k] - r2[k];
00050         meta_dist<k>::f(r_ij, r1, r2);
00051     }
00052 };
00053 
00054 template<>
00055 struct meta_dist<1> {
00056     static void f(double * const r_ij, double const * const r1,
00057             double const * const r2)
00058     {
00059         r_ij[0] = r1[0] - r2[0];
00060     }
00061 };  
00062 
00063 template<size_t ndim>
00064 struct cartesian_distance {
00065     static const size_t _ndim = ndim;
00066     void get_rij(double * const r_ij, double const * const r1,
00067             double const * const r2) const
00068     {
00069         static_assert(ndim > 0, "illegal box dimension");
00070         meta_dist<ndim>::f(r_ij, r1, r2);
00071     }
00072 };
00073 
00078 template<size_t IDX>
00079 struct meta_periodic_distance {
00080     static void f(double * const r_ij, double const * const r1,
00081                  double const * const r2, const double* _box, const double* _ibox)
00082     {
00083         const static size_t k = IDX - 1;
00084         r_ij[k] = r1[k] - r2[k];
00085         r_ij[k] -= round(r_ij[k] * _ibox[k]) * _box[k];
00086         meta_periodic_distance<k>::f(r_ij, r1, r2, _box, _ibox);
00087     } 
00088 };
00089 
00090 template<>
00091 struct meta_periodic_distance<1> {
00092     static void f(double * const r_ij, double const * const r1,
00093                  double const * const r2, const double* _box, const double* _ibox)
00094     {
00095         r_ij[0] = r1[0] - r2[0];
00096         r_ij[0] -= round(r_ij[0] * _ibox[0]) * _box[0];
00097     }
00098 };
00099 
00116 template<size_t IDX>
00117 struct meta_image {
00118     static void f(double *const x, const double* _ibox, const double* _box)
00119     {
00120         const static size_t k = IDX - 1;
00121         x[k] -= round(x[k] * _ibox[k]) * _box[k];
00122         meta_image<k>::f(x, _ibox, _box);
00123     }
00124 };
00125 
00126 template<>
00127 struct meta_image<1> {
00128     static void f(double *const x, const double* _ibox, const double* _box)
00129     {
00130         x[0] -= round(x[0] * _ibox[0]) * _box[0];
00131     }
00132 };
00133 
00134 template<size_t ndim>
00135 class periodic_distance {
00136 public:
00137     static const size_t _ndim = ndim;
00138     double _box[ndim];
00139     double _ibox[ndim];
00140 
00141     periodic_distance(Array<double> const box)
00142     {
00143         if (box.size() != _ndim) {
00144             throw std::invalid_argument("box.size() must be equal to ndim");
00145         }
00146         for (size_t i = 0; i < ndim; ++i) {
00147             _box[i] = box[i];
00148             _ibox[i] = 1 / box[i];
00149         }
00150     }
00151 
00152     periodic_distance()
00153     { 
00154         throw std::runtime_error("the empty constructor is not available for periodic boundaries"); 
00155     }
00156 
00157     inline void get_rij(double * const r_ij, double const * const r1,
00158                  double const * const r2) const
00159     {
00160         static_assert(ndim > 0, "illegal box dimension");
00161         meta_periodic_distance<ndim>::f(r_ij, r1, r2, _box, _ibox);
00162     }
00163 
00164     inline void put_in_box(Array<double>& coords) const
00165     {
00166         const size_t natoms = coords.size() / _ndim;
00167         for (size_t i = 0; i < natoms; ++i){
00168             const size_t i1 = i * _ndim;
00169             static_assert(ndim > 0, "illegal box dimension");
00170             meta_image<ndim>::f(&coords[i1], _ibox, _box);
00171         }
00172     }
00173 };
00174 
00175 /*
00176  * Interface classes to pass a generic non template pointer to DistanceInterface
00177  * (this should reduce the need for templating where best performance is not essential)
00178  */
00179 
00180 class DistanceInterface{
00181 protected:
00182 public:
00183     virtual void get_rij(double * const r_ij, double const * const r1, double const * const r2) const =0;
00184     virtual ~DistanceInterface(){ }
00185 };
00186 
00187 template<size_t ndim>
00188 class CartesianDistanceWrapper : public DistanceInterface{
00189 protected:
00190     cartesian_distance<ndim> _dist;
00191 public:
00192     static const size_t _ndim = ndim;
00193     inline void get_rij(double * const r_ij, double const * const r1,
00194             double const * const r2) const
00195     {
00196         _dist.get_rij(r_ij, r1, r2);
00197     }
00198 };
00199 
00200 template<size_t ndim>
00201 class PeriodicDistanceWrapper : public DistanceInterface{
00202 protected:
00203     periodic_distance<ndim> _dist;
00204 public:
00205     static const size_t _ndim = ndim;
00206     PeriodicDistanceWrapper(Array<double> const box)
00207         : _dist(box)
00208     {};
00209 
00210     inline void get_rij(double * const r_ij, double const * const r1,
00211             double const * const r2) const
00212     {
00213         _dist.get_rij(r_ij, r1, r2);
00214     }
00215 };
00216 
00217 } // namespace pele
00218 #endif // #ifndef _PELE_DISTANCE_H
 All Classes Namespaces Functions Variables Typedefs