Program Listing for File EvaluatorInversePower.hpp

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

// Hoomd based potentials
// Copyright (c) 2009-2022 The Regents of the University of Michigan.
// Part of HOOMD-blue, released under the BSD 3-Clause License.

#ifndef __PAIR_EVALUATOR_FORCE_SHIFTED_LJ_H__
#define __PAIR_EVALUATOR_FORCE_SHIFTED_LJ_H__

#ifndef __HIPCC__
#include <string>
#endif

#include "hoomd/HOOMDMath.h"
#include "hoomd/md/EvaluatorPairLJ.h"


// need to declare these class methods with __device__ qualifiers when building
// in nvcc DEVICE is __host__ __device__ when included in nvcc and blank when
// included into the host compiler
#ifdef __HIPCC__
#define DEVICE __device__
#define HOSTDEVICE __host__ __device__
#else
#define DEVICE
#define HOSTDEVICE
#endif

namespace hoomd_pele {
namespace md {

class EvaluatorPairForceShiftedLJ {
 public:
  typedef EvaluatorPairLJ::param_type param_type;

  DEVICE void load_shared(char *&ptr, unsigned int &available_bytes) {}

  HOSTDEVICE void allocate_shared(char *&ptr,
                                  unsigned int &available_bytes) const {}

#ifdef ENABLE_HIP
  void set_memory_hints() const {}
#endif


  DEVICE EvaluatorPairForceShiftedLJ(Scalar _rsq, Scalar _rcutsq,
                                     const param_type &_params)
      : rsq(_rsq),
        rcutsq(_rcutsq),
        lj1(_params.epsilon_x_4 * _params.sigma_6 * _params.sigma_6),
        lj2(_params.epsilon_x_4 * _params.sigma_6) {}

  DEVICE static bool needsDiameter() { return false; }

  DEVICE void setDiameter(Scalar di, Scalar dj) {}

  DEVICE static bool needsCharge() { return false; }

  DEVICE void setCharge(Scalar qi, Scalar qj) {}


  DEVICE bool evalForceAndEnergy(Scalar &force_divr, Scalar &pair_eng,
                                 bool energy_shift) {
    // compute the force divided by r in force_divr
    if (rsq < rcutsq && lj1 != 0) {
      Scalar r2inv = Scalar(1.0) / rsq;
      Scalar r6inv = r2inv * r2inv * r2inv;
      force_divr =
          r2inv * r6inv * (Scalar(12.0) * lj1 * r6inv - Scalar(6.0) * lj2);

      pair_eng = r6inv * (lj1 * r6inv - lj2);

      Scalar rcut2inv = Scalar(1.0) / rcutsq;
      Scalar rcut6inv = rcut2inv * rcut2inv * rcut2inv;

      if (energy_shift) pair_eng -= rcut6inv * (lj1 * rcut6inv - lj2);

      // shift force and add linear term to potential
      Scalar rcut_r_inv = fast::rsqrt(rsq * rcutsq);
      Scalar force_rcut_at_rcut =
          rcut6inv * (Scalar(12.0) * lj1 * rcut6inv - Scalar(6.0) * lj2);
      force_divr -= rcut_r_inv * force_rcut_at_rcut;
      pair_eng += (rsq * rcut_r_inv - Scalar(1.0)) * force_rcut_at_rcut;

      return true;
    } else
      return false;
  }

  DEVICE Scalar evalPressureLRCIntegral() { return 0; }

  DEVICE Scalar evalEnergyLRCIntegral() { return 0; }

#ifndef __HIPCC__

  static std::string getName() { return std::string("force_shift_lj"); }

  std::string getShapeSpec() const {
    throw std::runtime_error(
        "Shape definition not supported for this pair potential.");
  }
#endif

 protected:
  Scalar rsq;
  Scalar rcutsq;
  Scalar lj1;
  Scalar lj2;
};

}  // end namespace md
}  // end namespace hoomd_pele

#endif  // __PAIR_EVALUATOR_FORCE_SHIFTED_LJ_H__