libhoomd/computes/LJWallForceCompute.cc

   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 /*! \file LJWallForceCompute.cc
  53     \brief Defines the LJWallForceCompute class
  54 */
  55
  56 #ifdef WIN32
  57 #pragma warning( push )
  58 #pragma warning( disable : 4244 )
  59 #endif
  60
  61 #include <boost/python.hpp>
  62 using namespace boost::python;
  63
  64 #include "LJWallForceCompute.h"
  65 #include "WallData.h"
  66 #include <stdexcept>
  67
  68 using namespace std;
  69
  70 /*! \param sysdef System to compute forces on
  71     \param r_cut Cuttoff distance beyond which the force is zero.
  72 */
  73 LJWallForceCompute::LJWallForceCompute(boost::shared_ptr<SystemDefinition> sysdef, Scalar r_cut):
  74         ForceCompute(sysdef), m_r_cut(r_cut)
  75     {
  76     m_exec_conf->msg->notice(5) << "Constructing LJWallForceCompute" << endl;
  77
  78     if (r_cut < 0.0)
  79         {
  80         m_exec_conf->msg->error() << "wall.lj: Negative r_cut doesn't make sense." << endl;
  81         throw runtime_error("Error initializing LJWallForceCompute");
  82         }
  83
  84     // initialize the number of types value
  85     unsigned int ntypes = m_pdata->getNTypes();
  86     assert(ntypes > 0);
  87
  88     // allocate data for lj1 and lj2
  89     m_lj1 = new Scalar[ntypes];
  90     m_lj2 = new Scalar[ntypes];
  91
  92     assert(m_lj1);
  93     assert(m_lj2);
  94
  95     //initialize parameters to 0
  96     memset((void*)m_lj1, 0,sizeof(Scalar)*ntypes);
  97     memset((void*)m_lj2, 0,sizeof(Scalar)*ntypes);
  98
  99     }
 100
 101 /*! Frees used memory
 102 */
 103 LJWallForceCompute::~LJWallForceCompute()
 104     {
 105     m_exec_conf->msg->notice(5) << "Destroying LJWallForceCompute" << endl;
 106
 107     delete[] m_lj1;
 108     delete[] m_lj2;
 109     m_lj1 = NULL;
 110     m_lj2 = NULL;
 111     }
 112
 113 /*! \param typ Particle type index to set parameters for
 114     \param lj1 lj1 parameter
 115     \param lj2 lj2 parameter
 116
 117     \note \a lj1 are \a lj2 are low level parameters used in the calculation. In order to specify
 118     these for a normal lennard jones formula (with alpha), they should be set to the following.
 119     \a lj1 = 4.0 * epsilon * pow(sigma,12.0)
 120     \a lj2 = alpha * 4.0 * epsilon * pow(sigma,6.0);
 121 */
 122 void LJWallForceCompute::setParams(unsigned int typ, Scalar lj1, Scalar lj2)
 123     {
 124     if (typ >= m_pdata->getNTypes())
 125         {
 126         m_exec_conf->msg->error() << "wall.lj: Trying to set params for a non existant type! " << typ << endl;
 127         throw runtime_error("Error setting params in LJWallForceCompute");
 128         }
 129
 130     // set the parameters
 131     m_lj1[typ] = lj1;
 132     m_lj2[typ] = lj2;
 133     }
 134
 135 /*! LJWallForceCompute provides
 136     - \c wall_lj_energy
 137 */
 138 std::vector< std::string > LJWallForceCompute::getProvidedLogQuantities()
 139     {
 140     vector<string> list;
 141     list.push_back("wall_lj_energy");
 142     return list;
 143     }
 144
 145 Scalar LJWallForceCompute::getLogValue(const std::string& quantity, unsigned int timestep)
 146     {
 147     if (quantity == string("wall_lj_energy"))
 148         {
 149         compute(timestep);
 150         return calcEnergySum();
 151         }
 152     else
 153         {
 154         m_exec_conf->msg->error() << "wall.lj" << quantity << " is not a valid log quantity" << endl;
 155         throw runtime_error("Error getting log value");
 156         }
 157     }
 158
 159 void LJWallForceCompute::computeForces(unsigned int timestep)
 160     {
 161     // start the profile for this compute
 162     if (m_prof) m_prof->push("LJ wall");
 163
 164     // get numparticle var for easier access
 165     unsigned int numParticles = m_pdata->getN();
 166     boost::shared_ptr<WallData> wall_data = m_sysdef->getWallData();
 167     unsigned int numWalls = wall_data->getNumWalls();
 168
 169     // precalculate r_cut squqred
 170     Scalar r_cut_sq = m_r_cut * m_r_cut;
 171
 172     // precalculate box lengths for use in the periodic imaging
 173     BoxDim box = m_pdata->getBox();
 174
 175     // access the particle data
 176     ArrayHandle<Scalar4> h_pos(m_pdata->getPositions(), access_location::host, access_mode::read);
 177
 178     ArrayHandle<Scalar4> h_force(m_force,access_location::host, access_mode::overwrite);
 179     ArrayHandle<Scalar> h_virial(m_virial,access_location::host, access_mode::overwrite);
 180
 181     // Zero data for force calculation.
 182     memset((void*)h_force.data,0,sizeof(Scalar4)*m_force.getNumElements());
 183     memset((void*)h_virial.data,0,sizeof(Scalar)*m_virial.getNumElements());
 184
 185     // there are enough other checks on the input data: but it doesn't hurt to be safe
 186     assert(h_force.data);
 187     assert(h_virial.data);
 188
 189
 190     // here we go, main calc loop
 191     // loop over every particle in the sim,
 192     // calculate forces and store them in  m_force
 193     for (unsigned int i = 0; i < numParticles; i++)
 194         {
 195         // Initialize some force variables to be used as temporary
 196         // storage in each iteration, initialized to 0 from which the force will be computed
 197         Scalar3 f = make_scalar3(0, 0, 0);
 198         Scalar pe = 0.0;
 199
 200         // Grab particle data from all the arrays for this loop
 201         Scalar px = h_pos.data[i].x;
 202         Scalar py = h_pos.data[i].y;
 203         Scalar pz = h_pos.data[i].z;
 204         unsigned int type = __scalar_as_int(h_pos.data[i].w);
 205
 206         // for each wall that exists in the simulation
 207         // calculate the force that it exerts on a particle
 208         // the sum of the forces from each wall is the resulting force
 209         for (unsigned int cur_wall_idx = 0; cur_wall_idx < numWalls; cur_wall_idx++)
 210             {
 211             const Wall& cur_wall = wall_data->getWall(cur_wall_idx);
 212
 213             // calculate distance from point to plane
 214             // http://mathworld.wolfram.com/Point-PlaneDistance.html
 215             Scalar distFromWall = cur_wall.normal_x * (px - cur_wall.origin_x)
 216                                   + cur_wall.normal_y * (py - cur_wall.origin_y)
 217                                   + cur_wall.normal_z * (pz - cur_wall.origin_z);
 218
 219             // use the distance to create a vector pointing from the plane to the particle
 220             Scalar3 dx = make_scalar3(cur_wall.normal_x * distFromWall,
 221                                       cur_wall.normal_y * distFromWall,
 222                                       cur_wall.normal_z * distFromWall);
 223
 224             // continue with the evaluation of the LJ force copied from LJForceCompute
 225             // apply periodic boundary conditions
 226             dx = box.minImage(dx);
 227
 228             // start computing the force
 229             Scalar rsq = dot(dx,dx);
 230
 231             // only compute the force if the particles are closer than the cuttoff
 232             if (rsq < r_cut_sq)
 233                 {
 234                 // compute the force magnitude/r
 235                 Scalar r2inv = Scalar(1.0)/rsq;
 236                 Scalar r6inv = r2inv * r2inv * r2inv;
 237                 Scalar forcelj = r6inv * (Scalar(12.0)*m_lj1[type]*r6inv - Scalar(6.0)*m_lj2[type]);
 238                 Scalar fforce = forcelj * r2inv;
 239                 Scalar tmp_eng = r6inv * (m_lj1[type]*r6inv - m_lj2[type]);
 240
 241                 // accumulate the force vector
 242                 f += dx * fforce;
 243                 pe += tmp_eng;
 244                 }
 245             }
 246
 247         h_force.data[i].x = f.x;
 248         h_force.data[i].y = f.y;
 249         h_force.data[i].z = f.z;
 250         h_force.data[i].w = pe;
 251         }
 252
 253     if (m_prof) m_prof->pop();
 254     }
 255
 256 void export_LJWallForceCompute()
 257     {
 258     class_<LJWallForceCompute, boost::shared_ptr<LJWallForceCompute>, bases<ForceCompute>, boost::noncopyable >
 259     ("LJWallForceCompute", init< boost::shared_ptr<SystemDefinition>, Scalar >())
 260     .def("setParams", &LJWallForceCompute::setParams)
 261     ;
 262     }
 263
 264 #ifdef WIN32
 265 #pragma warning( pop )
 266 #endif