libhoomd/computes/CGCMMAngleForceCompute.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: dnlebard
  51
  52 #ifdef WIN32
  53 #pragma warning( push )
  54 #pragma warning( disable : 4244 )
  55 #endif
  56
  57 #include <boost/python.hpp>
  58 using namespace boost::python;
  59
  60 #include "CGCMMAngleForceCompute.h"
  61
  62 #include <iostream>
  63 #include <sstream>
  64 #include <stdexcept>
  65 #include <math.h>
  66
  67 using namespace std;
  68
  69 // \param SMALL a relatively small number
  70 #define SMALL Scalar(0.001)
  71
  72 /*! \file CGCMMAngleForceCompute.cc
  73     \brief Contains code for the CGCMMAngleForceCompute class
  74 */
  75
  76 /*! \param sysdef System to compute forces on
  77     \post Memory is allocated, and forces are zeroed.
  78 */
  79 CGCMMAngleForceCompute::CGCMMAngleForceCompute(boost::shared_ptr<SystemDefinition> sysdef)
  80     : ForceCompute(sysdef), m_K(NULL), m_t_0(NULL), m_eps(NULL), m_sigma(NULL), m_rcut(NULL), m_cg_type(NULL)
  81     {
  82     m_exec_conf->msg->notice(5) << "Constructing CGCMMAngleForceCompute" << endl;
  83
  84     // access the angle data for later use
  85     m_CGCMMAngle_data = m_sysdef->getAngleData();
  86
  87     // check for some silly errors a user could make
  88     if (m_CGCMMAngle_data->getNTypes() == 0)
  89         {
  90         m_exec_conf->msg->error() << "angle.cgcmm: No angle types specified" << endl;
  91         throw runtime_error("Error initializing CGCMMAngleForceCompute");
  92         }
  93
  94     // allocate the parameters
  95     m_K = new Scalar[m_CGCMMAngle_data->getNTypes()];
  96     m_t_0 = new Scalar[m_CGCMMAngle_data->getNTypes()];
  97     m_eps =  new Scalar[m_CGCMMAngle_data->getNTypes()];
  98     m_sigma = new Scalar[m_CGCMMAngle_data->getNTypes()];
  99     m_rcut =  new Scalar[m_CGCMMAngle_data->getNTypes()];
 100     m_cg_type = new unsigned int[m_CGCMMAngle_data->getNTypes()];
 101
 102     assert(m_K);
 103     assert(m_t_0);
 104     assert(m_eps);
 105     assert(m_sigma);
 106     assert(m_rcut);
 107     assert(m_cg_type);
 108
 109     memset((void*)m_K,0,sizeof(Scalar)*m_CGCMMAngle_data->getNTypes());
 110     memset((void*)m_t_0,0,sizeof(Scalar)*m_CGCMMAngle_data->getNTypes());
 111     memset((void*)m_eps,0,sizeof(Scalar)*m_CGCMMAngle_data->getNTypes());
 112     memset((void*)m_sigma,0,sizeof(Scalar)*m_CGCMMAngle_data->getNTypes());
 113     memset((void*)m_rcut,0,sizeof(Scalar)*m_CGCMMAngle_data->getNTypes());
 114     memset((void*)m_cg_type,0,sizeof(unsigned int)*m_CGCMMAngle_data->getNTypes());
 115
 116     prefact[0] = Scalar(0.0);
 117     prefact[1] = Scalar(6.75);
 118     prefact[2] = Scalar(2.59807621135332);
 119     prefact[3] = Scalar(4.0);
 120
 121     cgPow1[0]  = Scalar(0.0);
 122     cgPow1[1]  = Scalar(9.0);
 123     cgPow1[2]  = Scalar(12.0);
 124     cgPow1[3]  = Scalar(12.0);
 125
 126     cgPow2[0]  = Scalar(0.0);
 127     cgPow2[1]  = Scalar(6.0);
 128     cgPow2[2]  = Scalar(4.0);
 129     cgPow2[3]  = Scalar(6.0);
 130     }
 131
 132 CGCMMAngleForceCompute::~CGCMMAngleForceCompute()
 133     {
 134     m_exec_conf->msg->notice(5) << "Destroying CGCMMAngleForceCompute" << endl;
 135
 136     delete[] m_K;
 137     delete[] m_t_0;
 138     delete[] m_cg_type;
 139     delete[] m_eps;
 140     delete[] m_sigma;
 141     delete[] m_rcut;
 142     m_K = NULL;
 143     m_t_0 = NULL;
 144     m_cg_type = NULL;
 145     m_eps = NULL;
 146     m_sigma = NULL;
 147     m_rcut = NULL;
 148     }
 149
 150 /*! \param type Type of the angle to set parameters for
 151     \param K Stiffness parameter for the force computation
 152     \param t_0 Equilibrium angle in radians for the force computation
 153     \param cg_type the type of coarse grained angle
 154     \param eps the epsilon parameter for the 1-3 repulsion term
 155     \param sigma the sigma parameter for the 1-3 repulsion term
 156
 157     Sets parameters for the potential of a particular angle type
 158 */
 159 void CGCMMAngleForceCompute::setParams(unsigned int type, Scalar K, Scalar t_0, unsigned int cg_type, Scalar eps, Scalar sigma)
 160     {
 161     // make sure the type is valid
 162     if (type >= m_CGCMMAngle_data->getNTypes())
 163         {
 164         m_exec_conf->msg->error() << "angle.cgcmm: Invalid angle type specified" << endl;
 165         throw runtime_error("Error setting parameters in CGCMMAngleForceCompute");
 166         }
 167
 168     const double myPow1 = cgPow1[cg_type];
 169     const double myPow2 = cgPow2[cg_type];
 170
 171     Scalar my_rcut = sigma * ((Scalar) exp(1.0/(myPow1-myPow2)*log(myPow1/myPow2)));
 172
 173     m_K[type] = K;
 174     m_t_0[type] = t_0;
 175     m_cg_type[type] = cg_type;
 176     m_eps[type] = eps;
 177     m_sigma[type] = sigma;
 178     m_rcut[type] = my_rcut;
 179
 180     // check for some silly errors a user could make
 181     if (cg_type > 3)
 182         m_exec_conf->msg->warning() << "angle.cgcmm: Unrecognized exponents specified" << endl;
 183     if (K <= 0)
 184         m_exec_conf->msg->warning() << "angle.cgcmm: specified K <= 0" << endl;
 185     if (t_0 <= 0)
 186         m_exec_conf->msg->warning() << "angle.cgcmm: specified t_0 <= 0" << endl;
 187     if (eps <= 0)
 188         m_exec_conf->msg->warning() << "angle.cgcmm: specified eps <= 0" << endl;
 189     if (sigma <= 0)
 190         m_exec_conf->msg->warning() << "angle.cgcmm: specified sigma <= 0" << endl;
 191     }
 192
 193 /*! CGCMMAngleForceCompute provides
 194     - \c angle_cgcmm_energy
 195 */
 196 std::vector< std::string > CGCMMAngleForceCompute::getProvidedLogQuantities()
 197     {
 198     vector<string> list;
 199     list.push_back("angle_cgcmm_energy");
 200     return list;
 201     }
 202
 203 /*! \param quantity Name of the quantity to get the log value of
 204     \param timestep Current time step of the simulation
 205 */
 206 Scalar CGCMMAngleForceCompute::getLogValue(const std::string& quantity, unsigned int timestep)
 207     {
 208     if (quantity == string("angle_cgcmm_energy"))
 209         {
 210         compute(timestep);
 211         return calcEnergySum();
 212         }
 213     else
 214         {
 215         m_exec_conf->msg->error() << "angle.cgcmm: " << quantity << " is not a valid log quantity for CGCMMAngleForceCompute" << endl;
 216         throw runtime_error("Error getting log value");
 217         }
 218     }
 219
 220 /*! Actually perform the force computation
 221     \param timestep Current time step
 222  */
 223 void CGCMMAngleForceCompute::computeForces(unsigned int timestep)
 224     {
 225     if (m_prof) m_prof->push("CGCMMAngle");
 226
 227     assert(m_pdata);
 228     // access the particle data arrays
 229     ArrayHandle< unsigned int > h_rtag(m_pdata->getRTags(), access_location::host, access_mode::read);
 230     ArrayHandle< Scalar4 > h_pos(m_pdata->getPositions(), access_location::host, access_mode::read);
 231
 232     ArrayHandle<Scalar4> h_force(m_force,access_location::host, access_mode::overwrite);
 233     ArrayHandle<Scalar> h_virial(m_virial,access_location::host, access_mode::overwrite);
 234     unsigned int virial_pitch = m_virial.getPitch();
 235
 236     // Zero data for force calculation.
 237     memset((void*)h_force.data,0,sizeof(Scalar4)*m_force.getNumElements());
 238     memset((void*)h_virial.data,0,sizeof(Scalar)*m_virial.getNumElements());
 239
 240     // there are enough other checks on the input data: but it doesn't hurt to be safe
 241     assert(h_force.data);
 242     assert(h_virial.data);
 243     assert(h_pos.data);
 244     assert(h_rtag.data);
 245
 246     // get a local copy of the simulation box too
 247     const BoxDim& box = m_pdata->getGlobalBox();
 248
 249     // allocate forces
 250     Scalar fab[3], fcb[3];
 251     Scalar fac;
 252
 253     Scalar eac;
 254     Scalar vac[6];
 255     // for each of the angles
 256     const unsigned int size = (unsigned int)m_CGCMMAngle_data->getN();
 257     for (unsigned int i = 0; i < size; i++)
 258         {
 259         // lookup the tag of each of the particles participating in the angle
 260         const AngleData::members_t& angle = m_CGCMMAngle_data->getMembersByIndex(i);
 261         assert(angle.tag[0] < m_pdata->getNGlobal());
 262         assert(angle.tag[1] < m_pdata->getNGlobal());
 263         assert(angle.tag[1] < m_pdata->getNGlobal());
 264
 265         // transform a, b, and c into indicies into the particle data arrays
 266         // MEM TRANSFER: 6 ints
 267         unsigned int idx_a = h_rtag.data[angle.tag[0]];
 268         unsigned int idx_b = h_rtag.data[angle.tag[1]];
 269         unsigned int idx_c = h_rtag.data[angle.tag[2]];
 270
 271         // throw an error if this angle is incomplete
 272         if (idx_a == NOT_LOCAL|| idx_b == NOT_LOCAL || idx_c == NOT_LOCAL)
 273             {
 274             this->m_exec_conf->msg->error() << "angle.harmonic: angle " <<
 275                 angle.tag[0] << " " << angle.tag[1] << " " << angle.tag[2] << " incomplete." << endl << endl;
 276             throw std::runtime_error("Error in angle calculation");
 277             }
 278
 279         assert(idx_a < m_pdata->getN()+m_pdata->getNGhosts());
 280         assert(idx_b < m_pdata->getN()+m_pdata->getNGhosts());
 281         assert(idx_c < m_pdata->getN()+m_pdata->getNGhosts());
 282
 283         // calculate d\vec{r}
 284         // MEM_TRANSFER: 18 Scalars / FLOPS 9
 285         Scalar3 dab;
 286         dab.x = h_pos.data[idx_a].x - h_pos.data[idx_b].x;
 287         dab.y = h_pos.data[idx_a].y - h_pos.data[idx_b].y;
 288         dab.z = h_pos.data[idx_a].z - h_pos.data[idx_b].z;
 289
 290         Scalar3 dcb;
 291         dcb.x = h_pos.data[idx_c].x - h_pos.data[idx_b].x;
 292         dcb.y = h_pos.data[idx_c].y - h_pos.data[idx_b].y;
 293         dcb.z = h_pos.data[idx_c].z - h_pos.data[idx_b].z;
 294
 295         Scalar3 dac;
 296         dac.x = h_pos.data[idx_a].x - h_pos.data[idx_c].x; // used for the 1-3 JL interaction
 297         dac.y = h_pos.data[idx_a].y - h_pos.data[idx_c].y;
 298         dac.z = h_pos.data[idx_a].z - h_pos.data[idx_c].z;
 299
 300         // apply minimum image conventions to all 3 vectors
 301         dab = box.minImage(dab);
 302         dcb = box.minImage(dcb);
 303         dac = box.minImage(dac);
 304
 305         // on paper, the formula turns out to be: F = K*\vec{r} * (r_0/r - 1)
 306         // FLOPS: 14 / MEM TRANSFER: 2 Scalars
 307
 308         // FLOPS: 42 / MEM TRANSFER: 6 Scalars
 309         Scalar rsqab = dab.x*dab.x+dab.y*dab.y+dab.z*dab.z;
 310         Scalar rab = sqrt(rsqab);
 311         Scalar rsqcb = dcb.x*dcb.x+dcb.y*dcb.y+dcb.z*dcb.z;
 312         Scalar rcb = sqrt(rsqcb);
 313         Scalar rsqac = dac.x*dac.x+dac.y*dac.y+dac.z*dac.z;
 314         Scalar rac = sqrt(rsqac);
 315
 316         Scalar c_abbc = dab.x*dcb.x+dab.y*dcb.y+dab.z*dcb.z;
 317         c_abbc /= rab*rcb;
 318
 319         if (c_abbc > 1.0) c_abbc = 1.0;
 320         if (c_abbc < -1.0) c_abbc = -1.0;
 321
 322         Scalar s_abbc = sqrt(1.0 - c_abbc*c_abbc);
 323         if (s_abbc < SMALL) s_abbc = SMALL;
 324         s_abbc = 1.0/s_abbc;
 325
 326         //////////////////////////////////////////
 327         // THIS CODE DOES THE 1-3 LJ repulsions //
 328         //////////////////////////////////////////////////////////////////////////////
 329         fac = Scalar(0.0);
 330         eac = Scalar(0.0);
 331         for (int k = 0; k < 6; k++)
 332             vac[k] = Scalar(0.0);
 333
 334         unsigned int angle_type = m_CGCMMAngle_data->getTypeByIndex(i);
 335         if (rac < m_rcut[angle_type])
 336             {
 337             const unsigned int cg_type = m_cg_type[angle_type];
 338             const Scalar cg_pow1 = cgPow1[cg_type];
 339             const Scalar cg_pow2 = cgPow2[cg_type];
 340             const Scalar cg_pref = prefact[cg_type];
 341
 342             const Scalar cg_ratio = m_sigma[angle_type]/rac;
 343             const Scalar cg_eps   = m_eps[angle_type];
 344
 345             fac = cg_pref*cg_eps / rsqac * (cg_pow1*pow(cg_ratio,cg_pow1) - cg_pow2*pow(cg_ratio,cg_pow2));
 346             eac = cg_eps + cg_pref*cg_eps * (pow(cg_ratio,cg_pow1) - pow(cg_ratio,cg_pow2));
 347
 348             vac[0] = fac * dac.x*dac.x;
 349             vac[1] = fac * dac.x*dac.y;
 350             vac[2] = fac * dac.x*dac.z;
 351             vac[3] = fac * dac.y*dac.y;
 352             vac[4] = fac * dac.y*dac.z;
 353             vac[5] = fac * dac.z*dac.z;
 354             }
 355         //////////////////////////////////////////////////////////////////////////////
 356
 357         // actually calculate the force
 358         Scalar dth = acos(c_abbc) - m_t_0[angle_type];
 359         Scalar tk = m_K[angle_type]*dth;
 360
 361         Scalar a = -1.0 * tk * s_abbc;
 362         Scalar a11 = a*c_abbc/rsqab;
 363         Scalar a12 = -a / (rab*rcb);
 364         Scalar a22 = a*c_abbc / rsqcb;
 365
 366         fab[0] = a11*dab.x + a12*dcb.x;
 367         fab[1] = a11*dab.y + a12*dcb.y;
 368         fab[2] = a11*dab.z + a12*dcb.z;
 369
 370         fcb[0] = a22*dcb.x + a12*dab.x;
 371         fcb[1] = a22*dcb.y + a12*dab.y;
 372         fcb[2] = a22*dcb.z + a12*dab.z;
 373
 374         // compute 1/3 of the energy, 1/3 for each atom in the angle
 375         Scalar angle_eng = (0.5*tk*dth + eac)*Scalar(1.0/3.0);
 376
 377         // compute 1/3 of the virial, 1/3 for each atom in the angle
 378         // upper triangular version of virial tensor
 379         Scalar angle_virial[6];
 380         angle_virial[0] = Scalar(1./3.) * ( dab.x*fab[0] + dcb.x*fcb[0] );
 381         angle_virial[1] = Scalar(1./3.) * ( dab.y*fab[0] + dcb.y*fcb[0] );
 382         angle_virial[2] = Scalar(1./3.) * ( dab.z*fab[0] + dcb.z*fcb[0] );
 383         angle_virial[3] = Scalar(1./3.) * ( dab.y*fab[1] + dcb.y*fcb[1] );
 384         angle_virial[4] = Scalar(1./3.) * ( dab.z*fab[1] + dcb.z*fcb[1] );
 385         angle_virial[5] = Scalar(1./3.) * ( dab.z*fab[2] + dcb.z*fcb[2] );
 386         Scalar virial[6];
 387         for (unsigned int k=0; k < 6; k++)
 388             virial[k] = angle_virial[k] + Scalar(1./3.)*vac[k];
 389
 390         // Now, apply the force to each individual atom a,b,c, and accumlate the energy/virial
 391         // only apply force to local particles
 392         if (idx_a < m_pdata->getN())
 393             {
 394             h_force.data[idx_a].x += fab[0] + fac*dac.x;
 395             h_force.data[idx_a].y += fab[1] + fac*dac.y;
 396             h_force.data[idx_a].z += fab[2] + fac*dac.z;
 397             h_force.data[idx_a].w += angle_eng;
 398             for (int k = 0; k < 6; k++)
 399                 h_virial.data[k*virial_pitch+idx_a] += virial[k];
 400             }
 401
 402         if (idx_b < m_pdata->getN())
 403             {
 404             h_force.data[idx_b].x -= fab[0] + fcb[0];
 405             h_force.data[idx_b].y -= fab[1] + fcb[1];
 406             h_force.data[idx_b].z -= fab[2] + fcb[2];
 407             h_force.data[idx_b].w += angle_eng;
 408             for (int k = 0; k < 6; k++)
 409                 h_virial.data[k*virial_pitch+idx_b] += virial[k];
 410             }
 411
 412         if (idx_c < m_pdata->getN())
 413             {
 414             h_force.data[idx_c].x += fcb[0] - fac*dac.x;
 415             h_force.data[idx_c].y += fcb[1] - fac*dac.y;
 416             h_force.data[idx_c].z += fcb[2] - fac*dac.z;
 417             h_force.data[idx_c].w += angle_eng;
 418             for (int k = 0; k < 6; k++)
 419                 h_virial.data[k*virial_pitch+idx_c] += virial[k];
 420             }
 421         }
 422     if (m_prof) m_prof->pop();
 423     }
 424
 425 void export_CGCMMAngleForceCompute()
 426     {
 427     class_<CGCMMAngleForceCompute, boost::shared_ptr<CGCMMAngleForceCompute>, bases<ForceCompute>, boost::noncopyable >
 428     ("CGCMMAngleForceCompute", init< boost::shared_ptr<SystemDefinition> >())
 429     .def("setParams", &CGCMMAngleForceCompute::setParams)
 430     ;
 431     }
 432
 433 #ifdef WIN32
 434 #pragma warning( pop )
 435 #endif