share/hoomd/plugin_template_pair_ext/cppmodule/EvaluatorPairLJ2.h

   1 /*
   2 Highly Optimized Object-oriented Many-particle Dynamics -- Blue Edition
   3 (HOOMD-blue) Open Source Software License Copyright 2008-2011 Ames Laboratory
   4 Iowa State University and The Regents of the University of Michigan All rights
   5 reserved.
   6
   7 HOOMD-blue may contain modifications ("Contributions") provided, and to which
   8 copyright is held, by various Contributors who have granted The Regents of the
   9 University of Michigan the right to modify and/or distribute such Contributions.
  10
  11 You may redistribute, use, and create derivate works of HOOMD-blue, in source
  12 and binary forms, provided you abide by the following conditions:
  13
  14 * Redistributions of source code must retain the above copyright notice, this
  15 list of conditions, and the following disclaimer both in the code and
  16 prominently in any materials provided with the distribution.
  17
  18 * Redistributions in binary form must reproduce the above copyright notice, this
  19 list of conditions, and the following disclaimer in the documentation and/or
  20 other materials provided with the distribution.
  21
  22 * All publications and presentations based on HOOMD-blue, including any reports
  23 or published results obtained, in whole or in part, with HOOMD-blue, will
  24 acknowledge its use according to the terms posted at the time of submission on:
  25 http://codeblue.umich.edu/hoomd-blue/citations.html
  26
  27 * Any electronic documents citing HOOMD-Blue will link to the HOOMD-Blue website:
  28 http://codeblue.umich.edu/hoomd-blue/
  29
  30 * Apart from the above required attributions, neither the name of the copyright
  31 holder nor the names of HOOMD-blue's contributors may be used to endorse or
  32 promote products derived from this software without specific prior written
  33 permission.
  34
  35 Disclaimer
  36
  37 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS ``AS IS'' AND
  38 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  39 WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, AND/OR ANY
  40 WARRANTIES THAT THIS SOFTWARE IS FREE OF INFRINGEMENT ARE DISCLAIMED.
  41
  42 IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
  43 INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  44 BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  45 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  46 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
  47 OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
  48 ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  49 */
  50
  51 // Maintainer: joaander
  52
  53 #ifndef __PAIR_EVALUATOR_LJ_H__
  54 #define __PAIR_EVALUATOR_LJ_H__
  55
  56 #ifndef NVCC
  57 #include <string>
  58 #endif
  59
  60 // always need to include hoomd_config first
  61 #include <hoomd/hoomd_config.h>
  62 #include <hoomd/HOOMDMath.h>
  63
  64 /*! \file EvaluatorPairLJ2.h
  65     \brief Defines the pair evaluator class for LJ potentials
  66     \details As the prototypical example of a MD pair potential, this also serves as the primary documetnation and
  67     base reference for the implementation of pair evaluators.
  68 */
  69
  70 // need to declare these class methods with __device__ qualifiers when building in nvcc
  71 //! DEVICE is __host__ __device__ when included in nvcc and blank when included into the host compiler
  72 #ifdef NVCC
  73 #define DEVICE __device__
  74 #else
  75 #define DEVICE
  76 #endif
  77
  78 //! Class for evaluating the LJ pair potential
  79 /*! Note: This is identical to the pair potential EvaluatorPairLJ in hoomd. It is included here as an example.
  80     <b>General Overview</b>
  81
  82     EvaluatorPairLJ is a low level computation class that computes the LJ pair potential V(r). As the standard
  83     MD potential, it also serves as a well documented example of how to write additional pair potentials. "Standard"
  84     pair potentials in hoomd are all handled via the template class PotentialPair. PotentialPair takes a potential
  85     evaluator as a template argument. In this way, all the complicated data mangament and other details of computing
  86     the pair force and potential on every single particle is only written once in the template class and the difference
  87     V(r) potentials that can be calculated are simply handled with various evaluator classes. Template instantiation is
  88     equivalent to inlining code, so there is no performance loss.
  89
  90     In hoomd, a "standard" pair potential is defined as V(rsq, rcutsq, params, di, dj, qi, qj), where rsq is the squared
  91     distance between the two particles, rcutsq is the cuttoff radius at which the potential goes to 0, params is any
  92     number of per type-pair parameters, di, dj are the diameters of particles i and j, and qi, qj are the charges of
  93     particles i and j respectively.
  94
  95     Diameter and charge are not always needed by a given pair evaluator, so it must provide the functions
  96     needsDiameter() and needsCharge() which return boolean values signifying if they need those quantities or not. A
  97     false return value notifies PotentialPair that it need not even load those valuse from memory, boosting performance.
  98
  99     If needsDiameter() returns true, a setDiameter(Scalar di, Scalar dj) method will be called to set the two diameters.
 100     Similarly, if needsCharge() returns true, a setCharge(Scalar qi, Scalar qj) method will be called to set the two
 101     charges.
 102
 103     All other arguments are common among all pair potentials and passed into the constructor. Coefficients are handled
 104     in a special way: the pair evaluator class (and PotentialPair) manage only a single parameter variable for each
 105     type pair. Pair potentials that need more than 1 parameter can specify that their param_type be a compound
 106     structure and reference that. For coalesced read performance on G200 GPUs, it is highly recommended that param_type
 107     is one of the following types: Scalar, Scalar2, Scalar4.
 108
 109     The program flow will proceed like this: When a potential between a pair of particles is to be evaluated, a
 110     PairEvaluator is instantiated, passing the common parameters to the constructor and calling setDiameter() and/or
 111     setCharge() if need be. Then, the evalForceAndEnergy() method is called to evaluate the force and energy (more
 112     on that later). Thus, the evaluator must save all of the values it needs to compute the force and energy in member
 113     variables.
 114
 115     evalForceAndEnergy() makes the necessary computations and sets the out parameters with the computed values.
 116     Specifically after the method complets, \a force_divr must be set to the value
 117     \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
 118     \f$ V(r) - V(r_{\mathrm{cut}}) \f$ if \a energy_shift is true.
 119
 120     A pair potential evaluator class is also used on the GPU. So all of its members must be declared with the
 121     DEVICE keyword before them to mark them __device__ when compiling in nvcc and blank otherwise. If any other code
 122     needs to diverge between the host and device (i.e., to use a special math function like __powf on the device), it
 123     can similarly be put inside an ifdef NVCC block.
 124
 125     <b>LJ specifics</b>
 126
 127     EvaluatorPairLJ evaluates the function:
 128     \f[ V_{\mathrm{LJ}}(r) = 4 \varepsilon \left[ \left( \frac{\sigma}{r} \right)^{12} -
 129                                             \alpha \left( \frac{\sigma}{r} \right)^{6} \right] \f]
 130     broken up as follows for efficiency
 131     \f[ V_{\mathrm{LJ}}(r) = r^{-6} \cdot \left( 4 \varepsilon \sigma^{12} \cdot r^{-6} -
 132                                             4 \alpha \varepsilon \sigma^{6} \right) \f]
 133     . Similarly,
 134     \f[ -\frac{1}{r} \frac{\partial V_{\mathrm{LJ}}}{\partial r} = r^{-2} \cdot r^{-6} \cdot
 135             \left( 12 \cdot 4 \varepsilon \sigma^{12} \cdot r^{-6} - 6 \cdot 4 \alpha \varepsilon \sigma^{6} \right) \f]
 136
 137     The LJ potential does not need diameter or charge. Two parameters are specified and stored in a Scalar2. \a lj1 is
 138     placed in \a params.x and \a lj2 is in \a params.y.
 139
 140     These are related to the standard lj parameters sigma and epsilon by:
 141     - \a lj1 = 4.0 * epsilon * pow(sigma,12.0)
 142     - \a lj2 = alpha * 4.0 * epsilon * pow(sigma,6.0);
 143
 144 */
 145 class EvaluatorPairLJ2
 146     {
 147     public:
 148         //! Define the parameter type used by this pair potential evaluator
 149         typedef Scalar2 param_type;
 150
 151         //! Constructs the pair potential evaluator
 152         /*! \param _rsq Squared distance beteen the particles
 153             \param _rcutsq Sqauared distance at which the potential goes to 0
 154             \param _params Per type pair parameters of this potential
 155         */
 156         DEVICE EvaluatorPairLJ2(Scalar _rsq, Scalar _rcutsq, const param_type& _params)
 157             : rsq(_rsq), rcutsq(_rcutsq), lj1(_params.x), lj2(_params.y)
 158             {
 159             }
 160
 161         //! LJ doesn't use diameter
 162         DEVICE static bool needsDiameter() { return false; }
 163         //! Accept the optional diameter values
 164         /*! \param di Diameter of particle i
 165             \param dj Diameter of particle j
 166         */
 167         DEVICE void setDiameter(Scalar di, Scalar dj) { }
 168
 169         //! LJ doesn't use charge
 170         DEVICE static bool needsCharge() { return false; }
 171         //! Accept the optional diameter values
 172         /*! \param qi Charge of particle i
 173             \param qj Charge of particle j
 174         */
 175         DEVICE void setCharge(Scalar qi, Scalar qj) { }
 176
 177         //! Evaluate the force and energy
 178         /*! \param force_divr Output parameter to write the computed force divided by r.
 179             \param pair_eng Output parameter to write the computed pair energy
 180             \param energy_shift If true, the potential must be shifted so that V(r) is continuous at the cutoff
 181             \note There is no need to check if rsq < rcutsq in this method. Cutoff tests are performed
 182                   in PotentialPair.
 183
 184             \return True if they are evaluated or false if they are not because we are beyond the cuttoff
 185         */
 186         DEVICE bool evalForceAndEnergy(Scalar& force_divr, Scalar& pair_eng, bool energy_shift)
 187             {
 188             // compute the force divided by r in force_divr
 189             if (rsq < rcutsq && lj1 != 0)
 190                 {
 191                 Scalar r2inv = Scalar(1.0)/rsq;
 192                 Scalar r6inv = r2inv * r2inv * r2inv;
 193                 force_divr= r2inv * r6inv * (Scalar(12.0)*lj1*r6inv - Scalar(6.0)*lj2);
 194
 195                 pair_eng = r6inv * (lj1*r6inv - lj2);
 196
 197                 if (energy_shift)
 198                     {
 199                     Scalar rcut2inv = Scalar(1.0)/rcutsq;
 200                     Scalar rcut6inv = rcut2inv * rcut2inv * rcut2inv;
 201                     pair_eng -= rcut6inv * (lj1*rcut6inv - lj2);
 202                     }
 203                 return true;
 204                 }
 205             else
 206                 return false;
 207             }
 208
 209         #ifndef NVCC
 210         //! Get the name of this potential
 211         /*! \returns The potential name. Must be short and all lowercase, as this is the name energies will be logged as
 212             via analyze.log.
 213         */
 214         static std::string getName()
 215             {
 216             return std::string("lj");
 217             }
 218         #endif
 219
 220     protected:
 221         Scalar rsq;     //!< Stored rsq from the constructor
 222         Scalar rcutsq;  //!< Stored rcutsq from the constructor
 223         Scalar lj1;     //!< lj1 parameter extracted from the params passed to the constructor
 224         Scalar lj2;     //!< lj2 parameter extracted from the params passed to the constructor
 225     };
 226
 227
 228 #endif // __PAIR_EVALUATOR_LJ_H__
 229