diff --git a/polymd/EvaluatorPairPolydisperseLJ106.h b/polymd/EvaluatorPairPolydisperseLJ106.h new file mode 100644 index 0000000..3393c25 --- /dev/null +++ b/polymd/EvaluatorPairPolydisperseLJ106.h @@ -0,0 +1,198 @@ +// Copyright (c) 2009-2019 The Regents of the University of Michigan +// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. + + +// Maintainer: joaander + +#ifndef __PAIR_EVALUATOR_POLYDISPERSE_LJ106_H__ +#define __PAIR_EVALUATOR_POLYDISPERSE_LJ106_H__ + +#ifndef NVCC +#include +#endif + +#include "hoomd/HOOMDMath.h" + +/*! \file EvaluatorPairPolydisperseLJ106.h + \brief Defines the pair evaluator class for LJ106 potentials + \details As the prototypical example of a MD pair potential, this also serves as the primary documentation and + base reference for the implementation of pair evaluators. +*/ + +// 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 NVCC +#define DEVICE __device__ +#else +#define DEVICE +#endif + +//! Class for evaluating the LJ106 pair potential +/*! General Overview + + EvaluatorPairPolydisperseLJ106 is a low level computation class that computes the LJ106 pair potential V(r). As the standard + MD potential, it also serves as a well documented example of how to write additional pair potentials. "Standard" + pair potentials in hoomd are all handled via the template class PotentialPair. PotentialPair takes a potential + evaluator as a template argument. In this way, all the complicated data management and other details of computing + the pair force and potential on every single particle is only written once in the template class and the difference + V(r) potentials that can be calculated are simply handled with various evaluator classes. Template instantiation is + equivalent to inlining code, so there is no performance loss. + + In hoomd, a "standard" pair potential is defined as V(rsq, rcutsq, params, di, dj, qi, qj), where rsq is the squared + distance between the two particles, rcutsq is the cutoff radius at which the potential goes to 0, params is any + number of per type-pair parameters, di, dj are the diameters of particles i and j, and qi, qj are the charges of + particles i and j respectively. + + Diameter and charge are not always needed by a given pair evaluator, so it must provide the functions + needsDiameter() and needsCharge() which return boolean values signifying if they need those quantities or not. A + false return value notifies PotentialPair that it need not even load those values from memory, boosting performance. + + If needsDiameter() returns true, a setDiameter(Scalar di, Scalar dj) method will be called to set the two diameters. + Similarly, if needsCharge() returns true, a setCharge(Scalar qi, Scalar qj) method will be called to set the two + charges. + + All other arguments are common among all pair potentials and passed into the constructor. Coefficients are handled + in a special way: the pair evaluator class (and PotentialPair) manage only a single parameter variable for each + type pair. Pair potentials that need more than 1 parameter can specify that their param_type be a compound + structure and reference that. For coalesced read performance on G200 GPUs, it is highly recommended that param_type + is one of the following types: Scalar, Scalar2, Scalar4. + + The program flow will proceed like this: When a potential between a pair of particles is to be evaluated, a + PairEvaluator is instantiated, passing the common parameters to the constructor and calling setDiameter() and/or + setCharge() if need be. Then, the evalForceAndEnergy() method is called to evaluate the force and energy (more + on that later). Thus, the evaluator must save all of the values it needs to compute the force and energy in member + variables. + + evalForceAndEnergy() makes the necessary computations and sets the out parameters with the computed values. + Specifically after the method completes, \a force_divr must be set to the value + \f$ -\frac{1}{r}\frac{\partial V}{\partial r}\f$ and \a pair_eng must be set to the value \f$ V(r) \f$ if \a energy_shift is false or + \f$ V(r) - V(r_{\mathrm{cut}}) \f$ if \a energy_shift is true. + + A pair potential evaluator class is also used on the GPU. So all of its members must be declared with the + DEVICE keyword before them to mark them __device__ when compiling in nvcc and blank otherwise. If any other code + needs to diverge between the host and device (i.e., to use a special math function like __powf on the device), it + can similarly be put inside an ifdef NVCC block. + + LJ106 specifics + + EvaluatorPairPolydisperseLJ106 evaluates the function: + \f[ V_{\mathrm{LJ106}}(r) = 4 \varepsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - + \alpha \left( \frac{\sigma}{r} \right)^{6} \right] \f] + broken up as follows for efficiency + \f[ V_{\mathrm{LJ106}}(r) = r^{-6} \cdot \left( 4 \varepsilon \sigma^{12} \cdot r^{-6} - + 4 \alpha \varepsilon \sigma^{6} \right) \f] + . Similarly, + \f[ -\frac{1}{r} \frac{\partial V_{\mathrm{LJ106}}}{\partial r} = r^{-2} \cdot r^{-6} \cdot + \left( 12 \cdot 4 \varepsilon \sigma^{12} \cdot r^{-6} - 6 \cdot 4 \alpha \varepsilon \sigma^{6} \right) \f] + + The LJ106 potential does not need diameter or charge. Two parameters are specified and stored in a Scalar2. \a lj1 is + placed in \a params.x and \a lj2 is in \a params.y. + + These are related to the standard lj parameters sigma and epsilon by: + - \a lj1 = 4.0 * epsilon * pow(sigma,12.0) + - \a lj2 = alpha * 4.0 * epsilon * pow(sigma,6.0); + +*/ +class EvaluatorPairPolydisperseLJ106 + { + public: + //! Define the parameter type used by this pair potential evaluator + typedef Scalar3 param_type; + + //! Constructs the pair potential evaluator + /*! \param _rsq Squared distance between the particles + \param _rcutsq Squared distance at which the potential goes to 0 + \param _params Per type pair parameters of this potential + */ + DEVICE EvaluatorPairPolydisperseLJ106(Scalar _rsq, Scalar _rcutsq, const param_type& _params) + : rsq(_rsq), rcutsq(_rcutsq), v0(_params.x), eps(_params.y), scaledr_cut(_params.z), c0(Scalar(0.038758)), c1(Scalar(-0.0092432)), c2(Scalar(0.00058888)) + { + c0 = (Scalar(-21) + Scalar(10)*pow(scaledr_cut,4))*v0/pow(scaledr_cut,10); + c1 = (Scalar(35) - Scalar(15)*pow(scaledr_cut,4))*v0/pow(scaledr_cut,12); + c2 = (Scalar(-15) + Scalar(6)*pow(scaledr_cut,4))*v0/pow(scaledr_cut,14); + } + + //! LJ106 doesn't use diameter + DEVICE static bool needsDiameter() { return true; } + //! Accept the optional diameter values + /*! \param di Diameter of particle i + \param dj Diameter of particle j + */ + DEVICE void setDiameter(Scalar di, Scalar dj) + { + d_i = di; + d_j = dj; + } + + //! LJ106 doesn't use charge + DEVICE static bool needsCharge() { return false; } + //! Accept the optional diameter values + /*! \param qi Charge of particle i + \param qj Charge of particle j + */ + DEVICE void setCharge(Scalar qi, Scalar qj) { } + + //! Evaluate the force and energy + /*! \param force_divr Output parameter to write the computed force divided by r. + \param pair_eng Output parameter to write the computed pair energy + \param energy_shift If true, the potential must be shifted so that V(r) is continuous at the cutoff + \note There is no need to check if rsq < rcutsq in this method. Cutoff tests are performed + in PotentialPair. + + \return True if they are evaluated or false if they are not because we are beyond the cutoff + */ + DEVICE bool evalForceAndEnergy(Scalar& force_divr, Scalar& pair_eng, bool energy_shift) + { + Scalar sigma = 0.5*(d_i+d_j)*(1-eps*fabs(d_i-d_j)); + Scalar actualcutsq = scaledr_cut*scaledr_cut*sigma*sigma; + // compute the force divided by r in force_divr + if (rsq < actualcutsq && v0 != 0) + { + Scalar r2inv = sigma*sigma*Scalar(1.0)/rsq; + Scalar _rsq = Scalar(1.0)*rsq/(sigma*sigma); + Scalar r6inv = r2inv * r2inv * r2inv; + force_divr = (Scalar(10.0)*v0*r6inv*r6inv-Scalar(6.0)*v0*r2inv*r6inv-Scalar(2.0)*c1 -Scalar(4.0)*c2*_rsq)/(sigma*sigma); + + //No energy shift is needed + pair_eng = v0*(r2inv*r2inv*r2inv*r2inv*r2inv-r6inv)+c0+c1*_rsq+c2*_rsq*_rsq; + + return true; + } + else + return false; + } + + #ifndef NVCC + //! Get the name of this potential + /*! \returns The potential name. Must be short and all lowercase, as this is the name energies will be logged as + via analyze.log. + */ + static std::string getName() + { + return std::string("polydisperse_lj"); + } + std::string getShapeSpec() const + { + throw std::runtime_error("Shape definition not supported for this pair potential."); + } + + #endif + + protected: + Scalar rsq; //!< Stored rsq from the constructor + Scalar rcutsq; //!< Stored rcutsq from the constructor + Scalar d_i; //!< d_i diameter of particle i + Scalar d_j; //!< d_j diameter of particle j + + //Parameters to be read + Scalar v0; + Scalar eps; + Scalar scaledr_cut; + + //Additional parameters to be computed from the ones I read + Scalar c0; + Scalar c1; + Scalar c2; + }; + +#endif // __PAIR_EVALUATOR_POLYDISPERSE_LJ106_H__