pele
Python energy landscape explorer
/home/js850/projects/pele/source/pele/inversepower.h
00001 #ifndef _PELE_INVERSEPOWER_H
00002 #define _PELE_INVERSEPOWER_H
00003 
00004 #include <memory>
00005 
00006 #include "simple_pairwise_potential.h"
00007 #include "simple_pairwise_ilist.h"
00008 #include "atomlist_potential.h"
00009 #include "distance.h"
00010 #include "meta_pow.h"
00011 #include "cell_list_potential.h"
00012 
00013 namespace pele {
00014 
00040 struct InversePower_interaction {
00041     double const _eps;
00042     double const _pow;
00043     Array<double> const _radii;
00044 
00045     InversePower_interaction(double pow, double eps, Array<double> const radii)
00046         : _eps(eps),
00047           _pow(pow),
00048           _radii(radii.copy())
00049     {}
00050 
00051     /* calculate energy from distance squared */
00052     double energy(double r2, size_t atomi, size_t atomj) const
00053     {
00054         double E;
00055         double r0 = _radii[atomi] + _radii[atomj]; //sum of the hard core radii
00056         double r = std::sqrt(r2);
00057         if (r >= r0){
00058             E = 0.;
00059         } else {
00060             E = std::pow((1 -r/r0), _pow) * _eps/_pow;
00061         }
00062         return E;
00063     }
00064 
00065     /* calculate energy and gradient from distance squared, gradient is in -(dv/drij)/|rij| */
00066     double energy_gradient(double r2, double *gij, size_t atomi, size_t atomj) const
00067     {
00068         double E;
00069         double r0 = _radii[atomi] + _radii[atomj]; //sum of the hard core radii
00070         double r = std::sqrt(r2);
00071 
00072         if (r >= r0){
00073             E = 0.;
00074             *gij = 0;
00075         } else {
00076             double factor = std::pow((1 -r/r0), _pow) * _eps;
00077             E =  factor / _pow;
00078             *gij =  - factor / ((r-r0)*r);
00079         }
00080 
00081         return E;
00082     }
00083 
00084     double inline energy_gradient_hessian(double r2, double *gij, double *hij, size_t atomi, size_t atomj) const
00085     {
00086         double E;
00087         double r0 = _radii[atomi] + _radii[atomj]; //sum of the hard core radii
00088         double r = std::sqrt(r2);
00089         if (r >= r0){
00090             E = 0.;
00091             *gij = 0;
00092             *hij=0;
00093         } else {
00094             double factor = std::pow((1 -r/r0), _pow) * _eps;
00095             double denom = 1.0 / (r-r0);
00096             E =  factor / _pow;
00097             *gij =  - factor * denom / r ;
00098             *hij = (_pow-1) * factor * denom * denom;
00099         }
00100 
00101         return E;
00102     }
00103 };
00104 
00105 template<int POW>
00106 struct InverseIntPower_interaction {
00107     double const _eps;
00108     double const _pow;
00109     Array<double> const _radii;
00110 
00111     InverseIntPower_interaction(double eps, Array<double> const radii)
00112         : _eps(eps),
00113           _pow(POW),
00114           _radii(radii.copy())
00115     {}
00116 
00117     /* calculate energy from distance squared */
00118     double energy(double r2, size_t atomi, size_t atomj) const
00119     {
00120         double E;
00121         double r0 = _radii[atomi] + _radii[atomj]; //sum of the hard core radii
00122         double r = std::sqrt(r2);
00123         if (r >= r0){
00124             E = 0.;
00125         } else {
00126             //E = std::pow((1 -r/r0), _pow) * _eps/_pow;
00127             E = pos_int_pow<POW>(1 -r/r0) * _eps/_pow;
00128         }
00129         return E;
00130     }
00131 
00132     /* calculate energy and gradient from distance squared, gradient is in -(dv/drij)/|rij| */
00133     double energy_gradient(double r2, double *gij, size_t atomi, size_t atomj) const
00134     {
00135         double E;
00136         double r0 = _radii[atomi] + _radii[atomj]; //sum of the hard core radii
00137         double r = std::sqrt(r2);
00138 
00139         if (r >= r0){
00140             E = 0.;
00141             *gij = 0;
00142         } else {
00143             //double factor = std::pow((1 -r/r0), _pow) * _eps;
00144             double factor = pos_int_pow<POW>(1 -r/r0) * _eps;
00145             E =  factor / _pow;
00146             *gij =  - factor / ((r-r0)*r);
00147         }
00148 
00149         return E;
00150     }
00151 
00152     double inline energy_gradient_hessian(double r2, double *gij, double *hij, size_t atomi, size_t atomj) const
00153     {
00154         double E;
00155         double r0 = _radii[atomi] + _radii[atomj]; //sum of the hard core radii
00156         double r = std::sqrt(r2);
00157         if (r >= r0){
00158             E = 0.;
00159             *gij = 0;
00160             *hij=0;
00161         } else {
00162             //double factor = std::pow((1 -r/r0), _pow) * _eps;
00163             double factor = pos_int_pow<POW>(1 -r/r0) * _eps;
00164             double denom = 1.0 / (r-r0);
00165             E =  factor / _pow;
00166             *gij =  - factor * denom / r ;
00167             *hij = (_pow-1) * factor * denom * denom;
00168         }
00169 
00170         return E;
00171     }
00172 };
00173 
00174 template<int POW2>
00175 struct InverseHalfIntPower_interaction {
00176     double const _eps;
00177     double const _pow;
00178     Array<double> const _radii;
00179 
00180     InverseHalfIntPower_interaction(double eps, Array<double> const radii)
00181         : _eps(eps),
00182           _pow(0.5 * POW2),
00183           _radii(radii.copy())
00184     {}
00185 
00186     /* calculate energy from distance squared */
00187     double energy(double r2, size_t atomi, size_t atomj) const
00188     {
00189         double E;
00190         double r0 = _radii[atomi] + _radii[atomj]; //sum of the hard core radii
00191         double r = std::sqrt(r2);
00192         if (r >= r0){
00193             E = 0.;
00194         } else {
00195             //E = std::pow((1 -r/r0), _pow) * _eps/_pow;
00196             E = pos_half_int_pow<POW2>(1 -r/r0) * _eps/_pow;
00197         }
00198         return E;
00199     }
00200 
00201     /* calculate energy and gradient from distance squared, gradient is in -(dv/drij)/|rij| */
00202     double energy_gradient(double r2, double *gij, size_t atomi, size_t atomj) const
00203     {
00204         double E;
00205         double r0 = _radii[atomi] + _radii[atomj]; //sum of the hard core radii
00206         double r = std::sqrt(r2);
00207 
00208         if (r >= r0){
00209             E = 0.;
00210             *gij = 0;
00211         } else {
00212             //double factor = std::pow((1 -r/r0), _pow) * _eps;
00213             double factor = pos_half_int_pow<POW2>(1 -r/r0) * _eps;
00214             E =  factor / _pow;
00215             *gij =  - factor / ((r-r0)*r);
00216         }
00217 
00218         return E;
00219     }
00220 
00221     double inline energy_gradient_hessian(double r2, double *gij, double *hij, size_t atomi, size_t atomj) const
00222     {
00223         double E;
00224         double r0 = _radii[atomi] + _radii[atomj]; //sum of the hard core radii
00225         double r = std::sqrt(r2);
00226         if (r >= r0){
00227             E = 0.;
00228             *gij = 0;
00229             *hij=0;
00230         } else {
00231             //double factor = std::pow((1 -r/r0), _pow) * _eps;
00232             double factor = pos_half_int_pow<POW2>(1 -r/r0) * _eps;
00233             double denom = 1.0 / (r-r0);
00234             E =  factor / _pow;
00235             *gij =  - factor * denom / r ;
00236             *hij = (_pow-1) * factor * denom * denom;
00237         }
00238 
00239         return E;
00240     }
00241 };
00242 
00243 
00244 //
00245 // combine the components (interaction, looping method, distance function) into
00246 // defined classes
00247 //
00248 
00252 template <size_t ndim>
00253 class InversePower : public SimplePairwisePotential< InversePower_interaction, cartesian_distance<ndim> > {
00254 public:
00255     InversePower(double pow, double eps, pele::Array<double> const radii)
00256         : SimplePairwisePotential< InversePower_interaction, cartesian_distance<ndim> >(
00257                 std::make_shared<InversePower_interaction>(pow, eps, radii),
00258                 std::make_shared<cartesian_distance<ndim> >()
00259           )
00260     {}
00261 };
00262 
00263 template <size_t ndim>
00264 class InversePowerPeriodic : public SimplePairwisePotential< InversePower_interaction, periodic_distance<ndim> > {
00265 public:
00266     InversePowerPeriodic(double pow, double eps, pele::Array<double> const radii, pele::Array<double> const boxvec)
00267         : SimplePairwisePotential< InversePower_interaction, periodic_distance<ndim> >(
00268                 std::make_shared<InversePower_interaction>(pow, eps, radii),
00269                 std::make_shared<periodic_distance<ndim> >(boxvec)
00270           )
00271     {}
00272 };
00273 
00278 template <size_t ndim, int POW>
00279 class InverseIntPower : public SimplePairwisePotential< InverseIntPower_interaction<POW>, cartesian_distance<ndim> > {
00280 public:
00281     InverseIntPower(double eps, pele::Array<double> const radii)
00282         : SimplePairwisePotential< InverseIntPower_interaction<POW>, cartesian_distance<ndim> >(
00283                 std::make_shared<InverseIntPower_interaction<POW> >(eps, radii),
00284                 std::make_shared<cartesian_distance<ndim> >()
00285           )
00286     {}
00287 };
00288 
00289 template <size_t ndim, int POW>
00290 class InverseIntPowerPeriodic : public SimplePairwisePotential< InverseIntPower_interaction<POW>, periodic_distance<ndim> > {
00291 public:
00292     InverseIntPowerPeriodic(double eps, pele::Array<double> const radii, pele::Array<double> const boxvec)
00293         : SimplePairwisePotential< InverseIntPower_interaction<POW>, periodic_distance<ndim> >(
00294                 std::make_shared<InverseIntPower_interaction<POW> >(eps, radii),
00295                 std::make_shared<periodic_distance<ndim> >(boxvec)
00296           )
00297     {}
00298 };
00299 
00304 template <size_t ndim, int POW2>
00305 class InverseHalfIntPower : public SimplePairwisePotential< InverseHalfIntPower_interaction<POW2>, cartesian_distance<ndim> > {
00306 public:
00307     InverseHalfIntPower(double eps, pele::Array<double> const radii)
00308         : SimplePairwisePotential< InverseHalfIntPower_interaction<POW2>, cartesian_distance<ndim> >(
00309                 std::make_shared<InverseHalfIntPower_interaction<POW2> >(eps, radii),
00310                 std::make_shared<cartesian_distance<ndim> >()
00311           )
00312     {}
00313 };
00314 
00315 template <size_t ndim, int POW2>
00316 class InverseHalfIntPowerPeriodic : public SimplePairwisePotential< InverseHalfIntPower_interaction<POW2>, periodic_distance<ndim> > {
00317 public:
00318     InverseHalfIntPowerPeriodic(double eps, pele::Array<double> const radii, pele::Array<double> const boxvec)
00319         : SimplePairwisePotential< InverseHalfIntPower_interaction<POW2>, periodic_distance<ndim> >(
00320                 std::make_shared<InverseHalfIntPower_interaction<POW2> >(eps, radii),
00321                 std::make_shared<periodic_distance<ndim> >(boxvec)
00322           )
00323     {}
00324 };
00325 
00326 template <size_t ndim>
00327 class InversePowerCellLists : public CellListPotential< InversePower_interaction, cartesian_distance<ndim> > {
00328 public:
00329     InversePowerCellLists(double pow, double eps,
00330             pele::Array<double> const radii, pele::Array<double> const boxvec,
00331             const double rcut,
00332             const double ncellx_scale = 1.0)
00333         : CellListPotential< InversePower_interaction, cartesian_distance<ndim> >(
00334                 std::make_shared<InversePower_interaction>(pow, eps, radii),
00335                 std::make_shared<cartesian_distance<ndim> >(),
00336                 boxvec, rcut, ncellx_scale)
00337     {}
00338 };
00339 
00340 
00341 template <size_t ndim>
00342 class InversePowerPeriodicCellLists : public CellListPotential< InversePower_interaction, periodic_distance<ndim> > {
00343 public:
00344     InversePowerPeriodicCellLists(double pow, double eps,
00345             pele::Array<double> const radii, pele::Array<double> const boxvec,
00346             const double rcut,
00347             const double ncellx_scale = 1.0)
00348         : CellListPotential< InversePower_interaction, periodic_distance<ndim> >(
00349                 std::make_shared<InversePower_interaction>(pow, eps, radii),
00350                 std::make_shared<periodic_distance<ndim> >(boxvec),
00351                 boxvec, rcut, ncellx_scale)
00352     {}
00353 };
00354 
00355 } //namespace pele
00356 
00357 #endif //#ifndef _PELE_INVERSEPOWER_H
00358 
00359 
 All Classes Namespaces Functions Variables Typedefs