share/hoomd/plugin_template_evaluators_ext/cppmodule/EvaluatorBondHarmonicDPD.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 #ifndef __BOND_EVALUATOR_HARMONIC_DPD_H__
  51 #define __BOND_EVALUATOR_HARMONIC_DPD_H__
  52
  53 #ifndef NVCC
  54 #include <string>
  55 #endif
  56
  57 #include <hoomd/HOOMDMath.h>
  58
  59 /*! \file EvaluatorBondHarmonicDPD.h
  60     \brief Defines the bond evaluator class for composite harmonic bond
  61            + DPD pair potentials
  62 */
  63
  64 // need to declare these class methods with __device__ qualifiers when building in nvcc
  65 //! DEVICE is __host__ __device__ when included in nvcc and blank when included into the host compiler
  66 #ifdef NVCC
  67 #define DEVICE __device__
  68 #else
  69 #define DEVICE
  70 #endif
  71
  72 // call different optimized sqrt functions on the host / device
  73 //! RSQRT is rsqrtf when included in nvcc and 1.0 / sqrt(x) when included into the host compiler
  74 #ifdef NVCC
  75 #define RSQRT(x) rsqrtf( (x) )
  76 #else
  77 #define RSQRT(x) Scalar(1.0) / sqrt( (x) )
  78 #endif
  79
  80 //! Class for evaluating the harmonic bond potential
  81 class EvaluatorBondHarmonicDPD
  82     {
  83     public:
  84         //! Define the parameter type used by this bond potential evaluator
  85         typedef Scalar4 param_type;
  86
  87         //! Constructs the bond potential evaluator
  88         /*! \param _rsq Squared distance beteen the particles
  89             \param _params Per bond type parameters of this potential
  90         */
  91         DEVICE EvaluatorBondHarmonicDPD(Scalar _rsq, const param_type& _params)
  92             : rsq(_rsq),K(_params.x), r_0(_params.y), rcut(_params.z),
  93               a(_params.w)
  94             {
  95             }
  96
  97         //! Harmonic doesn't use diameter
  98         DEVICE static bool needsDiameter() { return false; }
  99
 100         //! Accept the optional diameter values
 101         /*! \param da Diameter of particle a
 102             \param db Diameter of particle b
 103         */
 104         DEVICE void setDiameter(Scalar da, Scalar db) { }
 105
 106         //! Harmonic doesn't use charge
 107         DEVICE static bool needsCharge() { return false; }
 108
 109         //! Accept the optional charge values
 110         /*! \param qa Charge of particle a
 111             \param qb Charge of particle b
 112         */
 113         DEVICE void setCharge(Scalar qa, Scalar qb) { }
 114
 115         //! Evaluate the force and energy
 116         /*! \param force_divr Output parameter to write the computed force divided by r.
 117             \param bond_eng Output parameter to write the computed bond energy
 118
 119             \return True if they are evaluated or false if the bond
 120                     energy is divergent
 121         */
 122         DEVICE bool evalForceAndEnergy(Scalar& force_divr, Scalar& bond_eng)
 123             {
 124             // evaluate harmonic bond potential
 125             Scalar r = sqrt(rsq);
 126             force_divr = K * (r_0 / r - Scalar(1.0));
 127             if (!isfinite(force_divr)) return false;
 128             bond_eng = Scalar(0.5) * K * (r_0 - r) * (r_0 - r);
 129
 130             if (rsq < rcut*rcut)
 131                 {
 132                 // evaluate DPD pair potential
 133                 Scalar rinv = RSQRT(rsq);
 134                 Scalar r = Scalar(1.0) / rinv;
 135                 Scalar rcutinv = Scalar(1.0)/ rcut;
 136
 137                 force_divr += a*(rinv - rcutinv);
 138                 bond_eng += a * (rcut - r) - Scalar(1.0/2.0) * a * rcutinv * (rcut*rcut - rsq);
 139                 }
 140
 141             return true;
 142             }
 143
 144         #ifndef NVCC
 145         //! Get the name of this potential
 146         /*! \returns The potential name. Must be short and all lowercase, as this is the name energies will be logged as
 147             via analyze.log.
 148         */
 149         static std::string getName()
 150             {
 151             return std::string("harmonic+dpd");
 152             }
 153         #endif
 154
 155     protected:
 156         Scalar rsq;        //!< Stored rsq from the constructor
 157         Scalar K;          //!< K parameter
 158         Scalar r_0;        //!< r_0 parameter
 159         Scalar rcut;       //!< Cutoff for pair potential
 160         Scalar a;           //!< a parameter for pair potential
 161     };
 162
 163
 164 #endif // __BOND_EVALUATOR_HARMONIC_DPD_H__