Program Listing for File lj.hpp

↰ Return to documentation for file (source/pele/lj.hpp)

/******************************+
 * This is an example on how to define a new c potential. The class
 * internally calls the lennard jones fortran function but manages all
 * the parameters.
 *
 * The python bindings for this potential are in lj.pyx
 *
 * For an alternative pure cython implementation for this interface
 * see _lj_cython.pyx
 *
 */

#ifndef PYGMIN_LJ_H
#define PYGMIN_LJ_H

#include <memory>

#include "atomlist_potential.hpp"
#include "base_interaction.hpp"
#include "distance.hpp"
#include "frozen_atoms.hpp"
#include "simple_pairwise_ilist.hpp"
#include "simple_pairwise_potential.hpp"

namespace pele {

struct lj_interaction : BaseInteraction {
 private:
  double const _C6, _C12;
  double const _6C6, _12C12;
  double const _42C6, _156C12;

 public:
  lj_interaction(double C6, double C12)
      : _C6(C6),
        _C12(C12),
        _6C6(6. * _C6),
        _12C12(12. * _C12),
        _42C6(42. * _C6),
        _156C12(156. * _C12) {}

  /* calculate energy from distance squared */
  double inline energy(double r2, const double dij) const {
    double ir2 = 1.0 / r2;
    double ir6 = ir2 * ir2 * ir2;
    double ir12 = ir6 * ir6;

    return -_C6 * ir6 + _C12 * ir12;
  }

  /* calculate energy and gradient from distance squared, gradient is in g/|rij|
   */
  double inline energy_gradient(double r2, double *gij,
                                const double dij) const {
    double ir2 = 1.0 / r2;
    double ir6 = ir2 * ir2 * ir2;
    double ir12 = ir6 * ir6;

    *gij = (_12C12 * ir12 - _6C6 * ir6) * ir2;
    return -_C6 * ir6 + _C12 * ir12;
  }

  double inline energy_gradient_hessian(double r2, double *gij, double *hij,
                                        const double dij) const {
    double ir2 = 1.0 / r2;
    double ir6 = ir2 * ir2 * ir2;
    double ir12 = ir6 * ir6;
    *gij = (_12C12 * ir12 - _6C6 * ir6) * ir2;
    *hij = (_156C12 * ir12 - _42C6 * ir6) * ir2;
    return -_C6 * ir6 + _C12 * ir12;
  }
};

//
// combine the components (interaction, looping method, distance function) into
// defined classes
//

class LJ : public SimplePairwisePotential<lj_interaction> {
 public:
  LJ(double C6, double C12)
      : SimplePairwisePotential<lj_interaction>(
            std::make_shared<lj_interaction>(C6, C12)) {}
};

class LJPeriodic
    : public SimplePairwisePotential<lj_interaction, periodic_distance<3>> {
 public:
  LJPeriodic(double C6, double C12, Array<double> const boxvec)
      : SimplePairwisePotential<lj_interaction, periodic_distance<3>>(
            std::make_shared<lj_interaction>(C6, C12),
            std::make_shared<periodic_distance<3>>(boxvec)) {}
};

class LJNeighborList : public SimplePairwiseNeighborList<lj_interaction> {
 public:
  LJNeighborList(Array<size_t> &ilist, double C6, double C12)
      : SimplePairwiseNeighborList<lj_interaction>(
            std::make_shared<lj_interaction>(C6, C12), ilist) {}
};
}  // namespace pele
#endif