libhoomd/potentials/EvaluatorPairSLJ.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 // Maintainer: joaander
  51
  52 #ifndef __PAIR_EVALUATOR_SLJ_H__
  53 #define __PAIR_EVALUATOR_SLJ_H__
  54
  55 #ifndef NVCC
  56 #include <string>
  57 #endif
  58
  59 #include "HOOMDMath.h"
  60
  61 /*! \file EvaluatorPairSLJ.h
  62     \brief Defines the pair evaluator class for shifted Lennard-Jones potentials
  63 */
  64
  65 // need to declare these class methods with __device__ qualifiers when building in nvcc
  66 // DEVICE is __host__ __device__ when included in nvcc and blank when included into the host compiler
  67 #ifdef NVCC
  68 #define DEVICE __device__
  69 #else
  70 #define DEVICE
  71 #endif
  72
  73 //! Class for evaluating the Gaussian pair potential
  74 /*! <b>General Overview</b>
  75
  76     See EvaluatorPairLJ
  77
  78     <b>SLJ specifics</b>
  79
  80     EvaluatorPairSLJ evaluates the function:
  81     \f{eqnarray*}
  82     V_{\mathrm{SLJ}}(r)  = & 4 \varepsilon \left[ \left( \frac{\sigma}{r - \Delta} \right)^{12} -
  83                            \left( \frac{\sigma}{r - \Delta} \right)^{6} \right] & r < (r_{\mathrm{cut}} + \Delta) \\
  84                          = & 0 & r \ge (r_{\mathrm{cut}} + \Delta) \\
  85     \f}
  86     where \f$ \Delta = (d_i + d_j)/2 - 1 \f$ and \f$ d_i \f$ is the diameter of particle \f$ i \f$.
  87
  88     The SLJ potential does not need charge, but does need diameter. Two parameters are specified and stored in a
  89     Scalar2. \a lj1 is placed in \a params.x and \a lj2 is in \a params.y.
  90
  91     These are related to the standard lj parameters sigma and epsilon by:
  92     - \a lj1 = 4.0 * epsilon * pow(sigma,12.0)
  93     - \a lj2 = 4.0 * epsilon * pow(sigma,6.0);
  94
  95     Due to the way that SLJ modifies the cutoff condition, it will not function properly with the xplor shifting mode.
  96 */
  97 class EvaluatorPairSLJ
  98     {
  99     public:
 100         //! Define the parameter type used by this pair potential evaluator
 101         typedef Scalar2 param_type;
 102
 103         //! Constructs the pair potential evaluator
 104         /*! \param _rsq Squared distance beteen the particles
 105             \param _rcutsq Sqauared distance at which the potential goes to 0
 106             \param _params Per type pair parameters of this potential
 107         */
 108         DEVICE EvaluatorPairSLJ(Scalar _rsq, Scalar _rcutsq, const param_type& _params)
 109             : rsq(_rsq), rcutsq(_rcutsq), lj1(_params.x), lj2(_params.y)
 110             {
 111             }
 112
 113         //! SLJ uses diameter
 114         DEVICE static bool needsDiameter() { return true; }
 115         //! Accept the optional diameter values
 116         /*! \param di Diameter of particle i
 117             \param dj Diameter of particle j
 118         */
 119         DEVICE void setDiameter(Scalar di, Scalar dj)
 120             {
 121             delta = (di + dj) / Scalar(2.0) - Scalar(1.0);
 122             }
 123
 124         //! SLJ doesn't use charge
 125         DEVICE static bool needsCharge() { return false; }
 126         //! Accept the optional diameter values
 127         /*! \param qi Charge of particle i
 128             \param qj Charge of particle j
 129         */
 130         DEVICE void setCharge(Scalar qi, Scalar qj) { }
 131
 132         //! Evaluate the force and energy
 133         /*! \param force_divr Output parameter to write the computed force divided by r.
 134             \param pair_eng Output parameter to write the computed pair energy
 135             \param energy_shift If true, the potential must be shifted so that V(r) is continuous at the cutoff
 136             \note There is no need to check if rsq < rcutsq in this method. Cutoff tests are performed
 137                   in PotentialPair.
 138
 139             \return True if they are evaluated or false if they are not because we are beyond the cuttoff
 140         */
 141         DEVICE bool evalForceAndEnergy(Scalar& force_divr, Scalar& pair_eng, bool energy_shift)
 142             {
 143             // precompute some quantities
 144             Scalar rinv = fast::rsqrt(rsq);
 145             Scalar r = Scalar(1.0) / rinv;
 146             Scalar rcutinv = fast::rsqrt(rcutsq);
 147             Scalar rcut = Scalar(1.0) / rcutinv;
 148
 149             // compute the force divided by r in force_divr
 150             if (r < (rcut + delta) && lj1 != 0)
 151                 {
 152                 Scalar rmd = r - delta;
 153                 Scalar rmdinv = Scalar(1.0) / rmd;
 154                 Scalar rmd2inv = rmdinv * rmdinv;
 155                 Scalar rmd6inv = rmd2inv * rmd2inv * rmd2inv;
 156                 force_divr= rinv * rmdinv * rmd6inv * (Scalar(12.0)*lj1*rmd6inv - Scalar(6.0)*lj2);
 157
 158                 pair_eng = rmd6inv * (lj1*rmd6inv - lj2);
 159
 160                 if (energy_shift)
 161                     {
 162                     Scalar rcut2inv = rcutinv * rcutinv;
 163                     Scalar rcut6inv = rcut2inv * rcut2inv * rcut2inv;
 164                     pair_eng -= rcut6inv * (lj1*rcut6inv - lj2);
 165                     }
 166                 return true;
 167                 }
 168             else
 169                 return false;
 170             }
 171
 172         #ifndef NVCC
 173         //! Get the name of this potential
 174         /*! \returns The potential name. Must be short and all lowercase, as this is the name energies will be logged as
 175             via analyze.log.
 176         */
 177         static std::string getName()
 178             {
 179             return std::string("slj");
 180             }
 181         #endif
 182
 183     protected:
 184         Scalar rsq;     //!< Stored rsq from the constructor
 185         Scalar rcutsq;  //!< Stored rcutsq from the constructor
 186         Scalar lj1;     //!< lj1 parameter extracted from the params passed to the constructor
 187         Scalar lj2;     //!< lj2 parameter extracted from the params passed to the constructor
 188         Scalar delta;   //!< Delta parameter extracted from the call to setDiameter
 189     };
 190
 191 #endif // __PAIR_EVALUATOR_SLJ_H__