share/hoomd/plugin_template_evaluators_ext/cppmodule/EvaluatorPairLJ2.h

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