libhoomd/computes/EvaluatorExternalPeriodic.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: jglaser
  51
  52 #ifndef __EVALUATOR_EXTERNAL_PERIODIC_H__
  53 #define __EVALUATOR_EXTERNAL_PERIODIC_H__
  54
  55 #ifndef NVCC
  56 #include <string>
  57 #endif
  58
  59 #include <math.h>
  60 #include "HOOMDMath.h"
  61 #include "BoxDim.h"
  62
  63 /*! \file EvaluatorExternalPeriodic.h
  64     \brief Defines the external potential evaluator to induce a periodic ordered phase
  65 */
  66
  67 // need to declare these class methods with __device__ qualifiers when building in nvcc
  68 // DEVICE is __host__ __device__ when included in nvcc and blank when included into the host compiler
  69 #ifdef NVCC
  70 #define DEVICE __device__
  71 #else
  72 #define DEVICE
  73 #endif
  74
  75 // SCALARASINT resolves to __scalar_as_int on the device and to __scalar_as_int on the host
  76 #ifdef NVCC
  77 #define SCALARASINT(x) __scalar_as_int(x)
  78 #else
  79 #define SCALARASINT(x) __scalar_as_int(x)
  80 #endif
  81
  82 //! Class for evaluating sphere constraints
  83 /*! <b>General Overview</b>
  84     EvaluatorExternalPeriodic is an evaluator to induce a periodic modulation on the concentration profile
  85     in the system, e.g. to generate a periodic phase in a system of diblock copolymers.
  86
  87     The external potential \f$V(\vec{r}) \f$ is implemented using the following formula:
  88
  89     \f[
  90     V(\vec{r}) = A * \tanh\left[\frac{1}{2 \pi p w} \cos\left(p \vec{b}_i\cdot\vec{r}\right)\right]
  91     \f]
  92
  93     where \f$A\f$ is the ordering parameter, \f$\vec{b}_i\f$ is the reciprocal lattice vector direction
  94     \f$i=0..2\f$, \f$p\f$ the periodicity and \f$w\f$ the interface width
  95     (relative to the distance \f$2\pi/|\mathbf{b_i}|\f$ between planes in the \f$i\f$-direction).
  96     The modulation is one-dimensional. It extends along the lattice vector \f$\mathbf{a}_i\f$ of the
  97     simulation cell.
  98 */
  99 class EvaluatorExternalPeriodic
 100     {
 101     public:
 102
 103         //! type of parameters this external potential accepts
 104         typedef Scalar4 param_type;
 105
 106         //! Constructs the constraint evaluator
 107         /*! \param X position of particle
 108             \param box box dimensions
 109             \param params per-type parameters of external potential
 110         */
 111         DEVICE EvaluatorExternalPeriodic(Scalar3 X, const BoxDim& box, const param_type& params)
 112             : m_pos(X),
 113               m_box(box)
 114             {
 115             m_index=  SCALARASINT(params.x);
 116             m_orderParameter = params.y;
 117             m_interfaceWidth = params.z;
 118             m_periodicity =SCALARASINT(params.w);
 119             }
 120
 121         //! Evaluate the force, energy and virial
 122         /*! \param F force vector
 123             \param energy value of the energy
 124             \param virial array of six scalars for the upper triangular virial tensor
 125         */
 126         DEVICE void evalForceEnergyAndVirial(Scalar3& F, Scalar& energy, Scalar* virial)
 127             {
 128             Scalar3 a2 = make_scalar3(0,0,0);
 129             Scalar3 a3 = make_scalar3(0,0,0);
 130
 131             F.x = Scalar(0.0);
 132             F.y = Scalar(0.0);
 133             F.z = Scalar(0.0);
 134             energy = Scalar(0.0);
 135
 136             // For this potential, since it uses scaled positions, the virial is always zero.
 137             for (unsigned int i = 0; i < 6; i++)
 138                 virial[i] = Scalar(0.0);
 139
 140             Scalar V_box = m_box.getVolume();
 141             // compute the vector pointing from P to V
 142             if (m_index == 0)
 143                 {
 144                 a2 = m_box.getLatticeVector(1);
 145                 a3 = m_box.getLatticeVector(2);
 146                 }
 147             else if (m_index == 1)
 148                 {
 149                 a2 = m_box.getLatticeVector(2);
 150                 a3 = m_box.getLatticeVector(0);
 151                 }
 152             else if (m_index == 2)
 153                 {
 154                 a2 = m_box.getLatticeVector(0);
 155                 a3 = m_box.getLatticeVector(1);
 156                 }
 157
 158             Scalar3 b = Scalar(2.0*M_PI)*make_scalar3(a2.y*a3.z-a2.z*a3.y,
 159                                                       a2.z*a3.x-a2.x*a3.z,
 160                                                       a2.x*a3.y-a2.y*a3.x)/V_box;
 161             Scalar clipParameter, arg, clipcos, tanH, sechSq;
 162
 163             Scalar3 q = b*m_periodicity;
 164             clipParameter   = Scalar(1.0)/Scalar(2.0*M_PI)/(m_periodicity*m_interfaceWidth);
 165             arg = dot(m_pos,q);
 166             clipcos = clipParameter*fast::cos(arg);
 167             tanH = tanhf(clipcos);
 168             sechSq = (Scalar(1.0) - tanH*tanH);
 169
 170             F = m_orderParameter*sechSq*clipParameter*fast::sin(arg)*q;
 171             energy = m_orderParameter*tanH;
 172             }
 173
 174         #ifndef NVCC
 175         //! Get the name of this potential
 176         /*! \returns The potential name. Must be short and all lowercase, as this is the name energies will be logged as
 177             via analyze.log.
 178         */
 179         static std::string getName()
 180             {
 181             return std::string("periodic");
 182             }
 183         #endif
 184
 185     protected:
 186         Scalar3 m_pos;                //!< particle position
 187         BoxDim m_box;                 //!< box dimensions
 188         unsigned int m_index;         //!< cartesian index of direction along which the lammellae should be orientied
 189         Scalar m_orderParameter;      //!< ordering parameter
 190         Scalar m_interfaceWidth;      //!< width of interface between lamellae (relative to box length)
 191         unsigned int m_periodicity;   //!< number of lamellae of each type
 192    };
 193
 194
 195 #endif // __EVALUATOR_EXTERNAL_LAMELLAR_H__