test/unit/test_cgcmm_force.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
  51 #ifdef WIN32
  52 #pragma warning( push )
  53 #pragma warning( disable : 4103 4244 )
  54 #endif
  55
  56 #include <iostream>
  57
  58 #include <boost/bind.hpp>
  59 #include <boost/function.hpp>
  60 #include <boost/shared_ptr.hpp>
  61
  62 #include "CGCMMForceCompute.h"
  63 #ifdef ENABLE_CUDA
  64 #include "CGCMMForceComputeGPU.h"
  65 #endif
  66
  67 #include "NeighborListBinned.h"
  68 #include "Initializers.h"
  69
  70 #include <math.h>
  71
  72 using namespace std;
  73 using namespace boost;
  74
  75 /*! \file cgcmm_force_test.cc
  76     \brief Implements unit tests for CGCMMForceCompute and descendants
  77     \ingroup unit_tests
  78 */
  79
  80 //! Name the unit test module
  81 #define BOOST_TEST_MODULE CGCMMForceTests
  82 #include "boost_utf_configure.h"
  83
  84 //! Typedef'd CGCMMForceCompute factory
  85 typedef boost::function<boost::shared_ptr<CGCMMForceCompute> (boost::shared_ptr<SystemDefinition> sysdef, boost::shared_ptr<NeighborList> nlist, Scalar r_cut)> cgcmmforce_creator;
  86
  87 //! Test the ability of the cgcmm LJ12-4 force compute to actually calucate forces
  88 void cgcmm_force_particle124_test(cgcmmforce_creator cgcmm_creator, boost::shared_ptr<ExecutionConfiguration> exec_conf)
  89     {
  90     // this 3-particle test subtly checks several conditions
  91     // the particles are arranged on the x axis,  1   2   3
  92     // such that 2 is inside the cuttoff radius of 1 and 3, but 1 and 3 are outside the cuttoff
  93     // of course, the buffer will be set on the neighborlist so that 3 is included in it
  94     // thus, this case tests the ability of the force summer to sum more than one force on
  95     // a particle and ignore a particle outside the radius
  96
  97     // periodic boundary conditions will be handeled in another test
  98     boost::shared_ptr<SystemDefinition> sysdef_3(new SystemDefinition(3, BoxDim(1000.0), 1, 0, 0, 0, 0, exec_conf));
  99     boost::shared_ptr<ParticleData> pdata_3 = sysdef_3->getParticleData();
 100
 101     pdata_3->setPosition(0,make_scalar3(0.0,0.0,0.0));
 102     pdata_3->setPosition(1,make_scalar3(pow(3.0,1.0/8.0),0.0,0.0));
 103     pdata_3->setPosition(2,make_scalar3(2.0*pow(3.0,1.0/8.0),0.0,0.0));
 104
 105     boost::shared_ptr<NeighborList> nlist_3(new NeighborList(sysdef_3, Scalar(1.3), Scalar(3.0)));
 106     boost::shared_ptr<CGCMMForceCompute> fc_3 = cgcmm_creator(sysdef_3, nlist_3, Scalar(1.3));
 107
 108     // first test: setup a sigma of 1.0 so that all forces will be 0
 109     Scalar epsilon = Scalar(1.15);
 110     Scalar sigma = Scalar(1.0);
 111     Scalar alpha = Scalar(1.0);
 112     Scalar lj1 = Scalar(2.598076) * epsilon * pow(sigma,Scalar(12.0));
 113     Scalar lj2 = Scalar(0.0);
 114     Scalar lj3 = Scalar(0.0);
 115     Scalar lj4 = -alpha * Scalar(2.598076) * epsilon * pow(sigma,Scalar(4.0));
 116     //Scalar lj1 = Scalar(0.0);
 117     //Scalar lj2 = alpha * Scalar(6.75) * epsilon * pow(sigma,Scalar(6.0));
 118     //Scalar lj3 = Scalar(6.75) * epsilon * pow(sigma,Scalar(9.0));
 119     //Scalar lj4 = Scalar(0.0);
 120
 121     fc_3->setParams(0,0,lj1,lj2,lj3,lj4);
 122
 123     // compute the forces
 124     fc_3->compute(0);
 125
 126     {
 127     GPUArray<Scalar4>& force_array_1 =  fc_3->getForceArray();
 128     GPUArray<Scalar>& virial_array_1 =  fc_3->getVirialArray();
 129     unsigned int pitch = fc_3->getVirialArray().getPitch();
 130     ArrayHandle<Scalar4> h_force_1(force_array_1,access_location::host,access_mode::read);
 131     ArrayHandle<Scalar> h_virial_1(virial_array_1,access_location::host,access_mode::read);
 132     MY_BOOST_CHECK_SMALL(h_force_1.data[0].x, tol);
 133     MY_BOOST_CHECK_SMALL(h_force_1.data[0].y, tol);
 134     MY_BOOST_CHECK_SMALL(h_force_1.data[0].z, tol);
 135     MY_BOOST_CHECK_CLOSE(h_force_1.data[0].w, -0.575, tol);
 136     MY_BOOST_CHECK_SMALL(h_virial_1.data[0*pitch+0]
 137                         +h_virial_1.data[3*pitch+0]
 138                         +h_virial_1.data[5*pitch+0], tol);
 139
 140     MY_BOOST_CHECK_SMALL(h_force_1.data[1].x, tol);
 141     MY_BOOST_CHECK_SMALL(h_force_1.data[1].y, tol);
 142     MY_BOOST_CHECK_SMALL(h_force_1.data[1].z, tol);
 143     MY_BOOST_CHECK_CLOSE(h_force_1.data[1].w, -1.15, tol);
 144     MY_BOOST_CHECK_SMALL(h_virial_1.data[0*pitch+1]
 145                         +h_virial_1.data[3*pitch+1]
 146                         +h_virial_1.data[5*pitch+1], tol);
 147
 148     MY_BOOST_CHECK_SMALL(h_force_1.data[2].x, tol);
 149     MY_BOOST_CHECK_SMALL(h_force_1.data[2].y, tol);
 150     MY_BOOST_CHECK_SMALL(h_force_1.data[2].z, tol);
 151     MY_BOOST_CHECK_CLOSE(h_force_1.data[2].w, -0.575, tol);
 152     MY_BOOST_CHECK_SMALL(h_virial_1.data[0*pitch+2]
 153                         +h_virial_1.data[3*pitch+2]
 154                         +h_virial_1.data[5*pitch+2], tol);
 155     }
 156
 157     // now change sigma and alpha so we can check that it is computing the right force
 158     sigma = Scalar(1.2); // < bigger sigma should push particle 0 left and particle 2 right
 159     alpha = Scalar(0.45);
 160     lj1 = Scalar(2.598076) * epsilon * pow(sigma,Scalar(12.0));
 161     lj2 = Scalar(0.0);
 162     lj3 = Scalar(0.0);
 163     lj4 = -alpha * Scalar(2.598076) * epsilon * pow(sigma,Scalar(4.0));
 164     fc_3->setParams(0,0,lj1,lj2,lj3,lj4);
 165     fc_3->compute(1);
 166
 167     {
 168     GPUArray<Scalar4>& force_array_2 =  fc_3->getForceArray();
 169     GPUArray<Scalar>& virial_array_2 =  fc_3->getVirialArray();
 170     ArrayHandle<Scalar4> h_force_2(force_array_2,access_location::host,access_mode::read);
 171     ArrayHandle<Scalar> h_virial_2(virial_array_2,access_location::host,access_mode::read);
 172     unsigned int pitch = virial_array_2.getPitch();
 173     MY_BOOST_CHECK_CLOSE(h_force_2.data[0].x, -48.0146523, tol);
 174     MY_BOOST_CHECK_SMALL(h_force_2.data[0].y, tol);
 175     MY_BOOST_CHECK_SMALL(h_force_2.data[0].z, tol);
 176     MY_BOOST_CHECK_CLOSE(h_force_2.data[0].w, 1.758563, tol);
 177     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_2.data[0*pitch+0]
 178                                  +h_virial_2.data[3*pitch+0]
 179                                  +h_virial_2.data[5*pitch+0]), 9.18042374, tol);
 180
 181     // center particle should still be a 0 force by symmetry
 182     MY_BOOST_CHECK_SMALL(h_force_2.data[1].x, tol);
 183     MY_BOOST_CHECK_SMALL(h_force_2.data[1].y, 1e-5);
 184     MY_BOOST_CHECK_SMALL(h_force_2.data[1].z, 1e-5);
 185     // there is still an energy and virial, though
 186     MY_BOOST_CHECK_CLOSE(h_force_2.data[1].w, 3.517125, tol);
 187     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_2.data[0*pitch+1]
 188                                  +h_virial_2.data[3*pitch+1]
 189                                  +h_virial_2.data[5*pitch+1]),18.3608475, tol);
 190
 191     MY_BOOST_CHECK_CLOSE(h_force_2.data[2].x, 48.0146561, tol);
 192     MY_BOOST_CHECK_SMALL(h_force_2.data[2].y, tol);
 193     MY_BOOST_CHECK_SMALL(h_force_2.data[2].z, tol);
 194     MY_BOOST_CHECK_CLOSE(h_force_2.data[2].w, 1.758563, tol);
 195     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_2.data[0*pitch+2]
 196                                  +h_virial_2.data[3*pitch+2]
 197                                  +h_virial_2.data[5*pitch+2]),9.18042374, tol);
 198     }
 199
 200     // swap the order of particles 0 ans 2 in memory to check that the force compute handles this properly
 201     {
 202     ArrayHandle<Scalar4> h_pos(pdata_3->getPositions(), access_location::host, access_mode::readwrite);
 203     ArrayHandle<unsigned int> h_tag(pdata_3->getTags(), access_location::host, access_mode::readwrite);
 204     ArrayHandle<unsigned int> h_rtag(pdata_3->getRTags(), access_location::host, access_mode::readwrite);
 205
 206     h_pos.data[2].x = h_pos.data[2].y = h_pos.data[2].z = 0.0;
 207     h_pos.data[0].x = Scalar(2.0*pow(3.0,1.0/8.0)); h_pos.data[0].y = h_pos.data[0].z = 0.0;
 208
 209     h_tag.data[0] = 2;
 210     h_tag.data[2] = 0;
 211     h_rtag.data[0] = 2;
 212     h_rtag.data[2] = 0;
 213     }
 214
 215     // notify the particle data that we changed the order
 216     pdata_3->notifyParticleSort();
 217
 218     // recompute the forces at the same timestep, they should be updated
 219     fc_3->compute(1);
 220
 221     {
 222     GPUArray<Scalar4>& force_array_3 =  fc_3->getForceArray();
 223     GPUArray<Scalar>& virial_array_3 =  fc_3->getVirialArray();
 224     ArrayHandle<Scalar4> h_force_3(force_array_3,access_location::host,access_mode::read);
 225     ArrayHandle<Scalar> h_virial_3(virial_array_3,access_location::host,access_mode::read);
 226     MY_BOOST_CHECK_CLOSE(h_force_3.data[0].x, 48.0146523, tol);
 227     MY_BOOST_CHECK_CLOSE(h_force_3.data[2].x, -48.0146523, tol);
 228     }
 229     }
 230
 231 //! Test the ability of the cgcmm LJ9-6 force compute to actually calucate forces
 232 void cgcmm_force_particle96_test(cgcmmforce_creator cgcmm_creator, boost::shared_ptr<ExecutionConfiguration> exec_conf)
 233     {
 234     // this 3-particle test subtly checks several conditions
 235     // the particles are arranged on the x axis,  1   2   3
 236     // such that 2 is inside the cuttoff radius of 1 and 3, but 1 and 3 are outside the cuttoff
 237     // of course, the buffer will be set on the neighborlist so that 3 is included in it
 238     // thus, this case tests the ability of the force summer to sum more than one force on
 239     // a particle and ignore a particle outside the radius
 240
 241     // periodic boundary conditions will be handeled in another test
 242     boost::shared_ptr<SystemDefinition> sysdef_3(new SystemDefinition(3, BoxDim(1000.0), 1, 0, 0, 0, 0, exec_conf));
 243     boost::shared_ptr<ParticleData> pdata_3 = sysdef_3->getParticleData();
 244
 245     {
 246     ArrayHandle<Scalar4> h_pos(pdata_3->getPositions(), access_location::host, access_mode::readwrite);
 247     h_pos.data[0].x = h_pos.data[0].y = h_pos.data[0].z = 0.0;
 248     h_pos.data[1].x = Scalar(pow(1.5,1.0/3.0)); h_pos.data[1].y = h_pos.data[1].z = 0.0;
 249     h_pos.data[2].x = Scalar(2.0*pow(1.5,1.0/3.0)); h_pos.data[2].y = h_pos.data[2].z = 0.0;
 250     }
 251     boost::shared_ptr<NeighborList> nlist_3(new NeighborList(sysdef_3, Scalar(1.3), Scalar(3.0)));
 252     boost::shared_ptr<CGCMMForceCompute> fc_3 = cgcmm_creator(sysdef_3, nlist_3, Scalar(1.3));
 253
 254     // first test: setup a sigma of 1.0 so that all forces will be 0
 255     Scalar epsilon = Scalar(1.15);
 256     Scalar sigma = Scalar(1.0);
 257     Scalar alpha = Scalar(1.0);
 258     Scalar lj1 = Scalar(0.0);
 259     Scalar lj2 = Scalar(6.75) * epsilon * pow(sigma,Scalar(9.0));
 260     Scalar lj3 = -alpha * Scalar(6.75) * epsilon * pow(sigma,Scalar(6.0));
 261     Scalar lj4 = Scalar(0.0);
 262
 263     fc_3->setParams(0,0,lj1,lj2,lj3,lj4);
 264
 265     // compute the forces
 266     fc_3->compute(0);
 267
 268     {
 269     GPUArray<Scalar4>& force_array_4 =  fc_3->getForceArray();
 270     GPUArray<Scalar>& virial_array_4 =  fc_3->getVirialArray();
 271     unsigned int pitch = virial_array_4.getPitch();
 272     ArrayHandle<Scalar4> h_force_4(force_array_4,access_location::host,access_mode::read);
 273     ArrayHandle<Scalar> h_virial_4(virial_array_4,access_location::host,access_mode::read);
 274     MY_BOOST_CHECK_SMALL(h_force_4.data[0].x, tol);
 275     MY_BOOST_CHECK_SMALL(h_force_4.data[0].y, tol);
 276     MY_BOOST_CHECK_SMALL(h_force_4.data[0].z, tol);
 277     MY_BOOST_CHECK_CLOSE(h_force_4.data[0].w, -0.575, tol);
 278     MY_BOOST_CHECK_SMALL(h_virial_4.data[0*pitch+0]
 279                         +h_virial_4.data[3*pitch+0]
 280                         +h_virial_4.data[5*pitch+0], tol);
 281
 282     MY_BOOST_CHECK_SMALL(h_force_4.data[1].x, tol);
 283     MY_BOOST_CHECK_SMALL(h_force_4.data[1].y, tol);
 284     MY_BOOST_CHECK_SMALL(h_force_4.data[1].z, tol);
 285     MY_BOOST_CHECK_CLOSE(h_force_4.data[1].w, -1.15, tol);
 286     MY_BOOST_CHECK_SMALL(h_virial_4.data[0*pitch+1]
 287                         +h_virial_4.data[3*pitch+1]
 288                         +h_virial_4.data[5*pitch+1], tol);
 289
 290     MY_BOOST_CHECK_SMALL(h_force_4.data[2].x, tol);
 291     MY_BOOST_CHECK_SMALL(h_force_4.data[2].y, tol);
 292     MY_BOOST_CHECK_SMALL(h_force_4.data[2].z, tol);
 293     MY_BOOST_CHECK_CLOSE(h_force_4.data[2].w, -0.575, tol);
 294     MY_BOOST_CHECK_SMALL(h_virial_4.data[0*pitch+2]
 295                         +h_virial_4.data[3*pitch+2]
 296                         +h_virial_4.data[5*pitch+2], tol);
 297
 298     }
 299
 300     // now change sigma and alpha so we can check that it is computing the right force
 301     sigma = Scalar(1.2); // < bigger sigma should push particle 0 left and particle 2 right
 302     alpha = Scalar(0.45);
 303     lj1 = Scalar(0.0);
 304     lj2 = Scalar(6.75) * epsilon * pow(sigma,Scalar(9.0));
 305     lj3 = -alpha * Scalar(6.75) * epsilon * pow(sigma,Scalar(6.0));
 306     lj4 = Scalar(0.0);
 307     fc_3->setParams(0,0,lj1,lj2,lj3,lj4);
 308     fc_3->compute(1);
 309
 310     {
 311     GPUArray<Scalar4>& force_array_5 =  fc_3->getForceArray();
 312     GPUArray<Scalar>& virial_array_5 =  fc_3->getVirialArray();
 313     unsigned int pitch = virial_array_5.getPitch();
 314     ArrayHandle<Scalar4> h_force_5(force_array_5,access_location::host,access_mode::read);
 315     ArrayHandle<Scalar> h_virial_5(virial_array_5,access_location::host,access_mode::read);
 316     MY_BOOST_CHECK_CLOSE(h_force_5.data[0].x, -69.00675, tol);
 317     MY_BOOST_CHECK_SMALL(h_force_5.data[0].y, tol);
 318     MY_BOOST_CHECK_SMALL(h_force_5.data[0].z, tol);
 319     MY_BOOST_CHECK_CLOSE(h_force_5.data[0].w, 3.615877, tol);
 320     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_5.data[0*pitch+0]
 321                         +h_virial_5.data[3*pitch+0]
 322                         +h_virial_5.data[5*pitch+0]), 13.1655083,tol);
 323
 324     // center particle should still be a 0 force by symmetry
 325     MY_BOOST_CHECK_SMALL(h_force_5.data[1].x, tol);
 326     MY_BOOST_CHECK_SMALL(h_force_5.data[1].y, 1e-5);
 327     MY_BOOST_CHECK_SMALL(h_force_5.data[1].z, 1e-5);
 328     // there is still an energy and virial, though
 329     MY_BOOST_CHECK_CLOSE(h_force_5.data[1].w, 7.231755, tol);
 330     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_5.data[0*pitch+1]
 331                         +h_virial_5.data[3*pitch+1]
 332                         +h_virial_5.data[5*pitch+1]), 26.3310165,tol);
 333
 334     MY_BOOST_CHECK_CLOSE(h_force_5.data[2].x, 69.00675, tol);
 335     MY_BOOST_CHECK_SMALL(h_force_5.data[2].y, tol);
 336     MY_BOOST_CHECK_SMALL(h_force_5.data[2].z, tol);
 337     MY_BOOST_CHECK_CLOSE(h_force_5.data[2].w, 3.615877, tol);
 338     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_5.data[0*pitch+2]
 339                         +h_virial_5.data[3*pitch+2]
 340                         +h_virial_5.data[5*pitch+2]), 13.1655083,tol);
 341     }
 342
 343     // swap the order of particles 0 ans 2 in memory to check that the force compute handles this properly
 344     {
 345     ArrayHandle<Scalar4> h_pos(pdata_3->getPositions(), access_location::host, access_mode::readwrite);
 346     ArrayHandle<unsigned int> h_tag(pdata_3->getTags(), access_location::host, access_mode::readwrite);
 347     ArrayHandle<unsigned int> h_rtag(pdata_3->getRTags(), access_location::host, access_mode::readwrite);
 348
 349     h_pos.data[2].x = h_pos.data[2].y = h_pos.data[2].z = 0.0;
 350     h_pos.data[0].x = Scalar(2.0*pow(1.5,1.0/3.0)); h_pos.data[0].y = h_pos.data[0].z = 0.0;
 351
 352     h_tag.data[0] = 2;
 353     h_tag.data[2] = 0;
 354     h_rtag.data[0] = 2;
 355     h_rtag.data[2] = 0;
 356     }
 357
 358     // notify the particle data that we changed the order
 359     pdata_3->notifyParticleSort();
 360
 361     // recompute the forces at the same timestep, they should be updated
 362     fc_3->compute(1);
 363
 364     {
 365     GPUArray<Scalar4>& force_array_6 =  fc_3->getForceArray();
 366     GPUArray<Scalar>& virial_array_6 =  fc_3->getVirialArray();
 367     ArrayHandle<Scalar4> h_force_6(force_array_6,access_location::host,access_mode::read);
 368     ArrayHandle<Scalar> h_virial_6(virial_array_6,access_location::host,access_mode::read);
 369     MY_BOOST_CHECK_CLOSE(h_force_6.data[0].x, 69.00675, tol);
 370     MY_BOOST_CHECK_CLOSE(h_force_6.data[2].x, -69.00675, tol);
 371     }
 372     }
 373
 374 //! Tests the ability of a CGCMMForceCompute to handle periodic boundary conditions
 375 void cgcmm_force_periodic_test(cgcmmforce_creator cgcmm_creator, boost::shared_ptr<ExecutionConfiguration> exec_conf)
 376     {
 377     ////////////////////////////////////////////////////////////////////
 378     // now, lets do a more thorough test and include boundary conditions
 379     // there are way too many permutations to test here, so I will simply
 380     // test +x, -x, +y, -y, +z, and -z independantly
 381     // build a 6 particle system with particles across each boundary
 382     // also test the ability of the force compute to use different particle types
 383
 384     boost::shared_ptr<SystemDefinition> sysdef_6(new SystemDefinition(6, BoxDim(20.0, 40.0, 60.0), 3, 0, 0, 0, 0, exec_conf));
 385     boost::shared_ptr<ParticleData> pdata_6 = sysdef_6->getParticleData();
 386
 387     pdata_6->setPosition(0,make_scalar3(-9.6,0.0,0.0));
 388     pdata_6->setPosition(1,make_scalar3(9.6,0.0,0.0));
 389     pdata_6->setPosition(2,make_scalar3(0.0,-19.6,0.0));
 390     pdata_6->setPosition(3,make_scalar3(0.0,19.6,0.0));
 391     pdata_6->setPosition(4,make_scalar3(0.0,0.0,-29.6));
 392     pdata_6->setPosition(5,make_scalar3(0.0,0.0,29.6));
 393
 394     pdata_6->setType(0,0);
 395     pdata_6->setType(1,1);
 396     pdata_6->setType(2,2);
 397     pdata_6->setType(3,0);
 398     pdata_6->setType(4,2);
 399     pdata_6->setType(5,1);
 400
 401     boost::shared_ptr<NeighborList> nlist_6(new NeighborList(sysdef_6, Scalar(1.3), Scalar(3.0)));
 402     boost::shared_ptr<CGCMMForceCompute> fc_6 = cgcmm_creator(sysdef_6, nlist_6, Scalar(1.3));
 403
 404     // choose a small sigma so that all interactions are attractive
 405     Scalar epsilon = Scalar(1.0);
 406     Scalar sigma = Scalar(0.5);
 407     Scalar alpha = Scalar(0.45);
 408     Scalar lj1 = Scalar(4.0) * epsilon * pow(sigma,Scalar(12.0));
 409     Scalar lj2 = Scalar(0.0);
 410     Scalar lj3 = -alpha * Scalar(4.0) * epsilon * pow(sigma,Scalar(6.0));
 411     Scalar lj4 = Scalar(0.0);
 412
 413     // make life easy: just change epsilon for the different pairs
 414     fc_6->setParams(0,0,lj1,lj2,lj3,lj4);
 415     fc_6->setParams(0,1,Scalar(2.0)*lj1,Scalar(2.0)*lj2,Scalar(2.0)*lj3,Scalar(2.0)*lj4);
 416     fc_6->setParams(0,2,Scalar(3.0)*lj1,Scalar(3.0)*lj2,Scalar(3.0)*lj3,Scalar(3.0)*lj4);
 417     fc_6->setParams(1,1,Scalar(4.0)*lj1,Scalar(4.0)*lj2,Scalar(4.0)*lj3,Scalar(4.0)*lj4);
 418     fc_6->setParams(1,2,Scalar(5.0)*lj1,Scalar(5.0)*lj2,Scalar(5.0)*lj3,Scalar(5.0)*lj4);
 419     fc_6->setParams(2,2,Scalar(6.0)*lj1,Scalar(6.0)*lj2,Scalar(6.0)*lj3,Scalar(6.0)*lj4);
 420
 421     fc_6->compute(0);
 422
 423     {
 424     GPUArray<Scalar4>& force_array_7 =  fc_6->getForceArray();
 425     GPUArray<Scalar>& virial_array_7 =  fc_6->getVirialArray();
 426     unsigned int pitch = virial_array_7.getPitch();
 427     ArrayHandle<Scalar4> h_force_7(force_array_7,access_location::host,access_mode::read);
 428     ArrayHandle<Scalar> h_virial_7(virial_array_7,access_location::host,access_mode::read);
 429     // particle 0 should be pulled left
 430     MY_BOOST_CHECK_CLOSE(h_force_7.data[0].x, -1.18299976747949, tol);
 431     MY_BOOST_CHECK_SMALL(h_force_7.data[0].y, tol);
 432     MY_BOOST_CHECK_SMALL(h_force_7.data[0].z, tol);
 433     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_7.data[0*pitch+0]
 434                                  +h_virial_7.data[3*pitch+0]
 435                                  +h_virial_7.data[5*pitch+0]), -0.15773330233059, tol);
 436
 437     // particle 1 should be pulled right
 438     MY_BOOST_CHECK_CLOSE(h_force_7.data[1].x, 1.18299976747949, tol);
 439     MY_BOOST_CHECK_SMALL(h_force_7.data[1].y, 1e-5);
 440     MY_BOOST_CHECK_SMALL(h_force_7.data[1].z, 1e-5);
 441     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_7.data[0*pitch+1]
 442                                  +h_virial_7.data[3*pitch+1]
 443                                  +h_virial_7.data[5*pitch+1]), -0.15773330233059, tol);
 444
 445     // particle 2 should be pulled down
 446     MY_BOOST_CHECK_CLOSE(h_force_7.data[2].y, -1.77449965121923, tol);
 447     MY_BOOST_CHECK_SMALL(h_force_7.data[2].x, tol);
 448     MY_BOOST_CHECK_SMALL(h_force_7.data[2].z, tol);
 449     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_7.data[0*pitch+2]
 450                                  +h_virial_7.data[3*pitch+2]
 451                                  +h_virial_7.data[5*pitch+2]), -0.23659995349591, tol);
 452
 453     // particle 3 should be pulled up
 454     MY_BOOST_CHECK_CLOSE(h_force_7.data[3].y, 1.77449965121923, tol);
 455     MY_BOOST_CHECK_SMALL(h_force_7.data[3].x, 1e-5);
 456     MY_BOOST_CHECK_SMALL(h_force_7.data[3].z, 1e-5);
 457     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_7.data[0*pitch+3]
 458                                  +h_virial_7.data[3*pitch+3]
 459                                  +h_virial_7.data[5*pitch+3]), -0.23659995349591, tol);
 460
 461     // particle 4 should be pulled back
 462     MY_BOOST_CHECK_CLOSE(h_force_7.data[4].z, -2.95749941869871, tol);
 463     MY_BOOST_CHECK_SMALL(h_force_7.data[4].x, tol);
 464     MY_BOOST_CHECK_SMALL(h_force_7.data[4].y, tol);
 465     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_7.data[0*pitch+4]
 466                                  +h_virial_7.data[3*pitch+4]
 467                                  +h_virial_7.data[5*pitch+4]),  -0.39433325582651, tol);
 468
 469     // particle 3 should be pulled forward
 470     MY_BOOST_CHECK_CLOSE(h_force_7.data[5].z, 2.95749941869871, tol);
 471     MY_BOOST_CHECK_SMALL(h_force_7.data[5].x, 1e-5);
 472     MY_BOOST_CHECK_SMALL(h_force_7.data[5].y, 1e-5);
 473     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_7.data[0*pitch+5]
 474                                  +h_virial_7.data[3*pitch+5]
 475                                  +h_virial_7.data[5*pitch+5]),  -0.39433325582651, tol);
 476     }
 477     }
 478
 479 //! Unit test a comparison between 2 CGCMMForceComputes on a "real" system
 480 void cgcmm_force_comparison_test(cgcmmforce_creator cgcmm_creator1, cgcmmforce_creator cgcmm_creator2, boost::shared_ptr<ExecutionConfiguration> exec_conf)
 481     {
 482     const unsigned int N = 5000;
 483
 484     // create a random particle system to sum forces on
 485     RandomInitializer rand_init(N, Scalar(0.2), Scalar(0.9), "A");
 486     boost::shared_ptr<SnapshotSystemData> snap;
 487     snap = rand_init.getSnapshot();
 488     boost::shared_ptr<SystemDefinition> sysdef(new SystemDefinition(snap, exec_conf));
 489     boost::shared_ptr<NeighborListBinned> nlist(new NeighborListBinned(sysdef, Scalar(3.0), Scalar(0.8)));
 490
 491     boost::shared_ptr<CGCMMForceCompute> fc1 = cgcmm_creator1(sysdef, nlist, Scalar(3.0));
 492     boost::shared_ptr<CGCMMForceCompute> fc2 = cgcmm_creator2(sysdef, nlist, Scalar(3.0));
 493
 494     // setup some values for alpha and sigma
 495     Scalar epsilon = Scalar(1.0);
 496     Scalar sigma = Scalar(1.2);
 497     Scalar alpha = Scalar(0.45);
 498     Scalar lj1 = Scalar(4.0) * epsilon * pow(sigma,Scalar(12.0));
 499     Scalar lj2 = Scalar(0.0);
 500     Scalar lj3 = -alpha * Scalar(4.0) * epsilon * pow(sigma,Scalar(6.0));
 501     Scalar lj4 = Scalar(0.0);
 502
 503     // specify the force parameters
 504     fc1->setParams(0,0,lj1,lj2,lj3,lj4);
 505     fc2->setParams(0,0,lj1,lj2,lj3,lj4);
 506
 507     // compute the forces
 508     fc1->compute(0);
 509     fc2->compute(0);
 510
 511     {
 512     // verify that the forces are identical (within roundoff errors)
 513     GPUArray<Scalar4>& force_array_8 =  fc1->getForceArray();
 514     GPUArray<Scalar>& virial_array_8 =  fc1->getVirialArray();
 515     unsigned int pitch = virial_array_8.getPitch();
 516     ArrayHandle<Scalar4> h_force_8(force_array_8,access_location::host,access_mode::read);
 517     ArrayHandle<Scalar> h_virial_8(virial_array_8,access_location::host,access_mode::read);
 518     GPUArray<Scalar4>& force_array_9 =  fc2->getForceArray();
 519     GPUArray<Scalar>& virial_array_9 =  fc2->getVirialArray();
 520     ArrayHandle<Scalar4> h_force_9(force_array_9,access_location::host,access_mode::read);
 521     ArrayHandle<Scalar> h_virial_9(virial_array_9,access_location::host,access_mode::read);
 522
 523     // compare average deviation between the two computes
 524     double deltaf2 = 0.0;
 525     double deltape2 = 0.0;
 526     double deltav2[6];
 527     for (unsigned int i = 0; i < 6; i++)
 528         deltav2[i] = 0.0;
 529
 530     for (unsigned int i = 0; i < N; i++)
 531         {
 532         deltaf2 += double(h_force_9.data[i].x - h_force_8.data[i].x) * double(h_force_9.data[i].x - h_force_8.data[i].x);
 533         deltaf2 += double(h_force_9.data[i].y - h_force_8.data[i].y) * double(h_force_9.data[i].y - h_force_8.data[i].y);
 534         deltaf2 += double(h_force_9.data[i].z - h_force_8.data[i].z) * double(h_force_9.data[i].z - h_force_8.data[i].z);
 535         deltape2 += double(h_force_9.data[i].w - h_force_8.data[i].w) * double(h_force_9.data[i].w - h_force_8.data[i].w);
 536         for (int j = 0; j < 6; j++)
 537             deltav2[j] += double(h_virial_9.data[j*pitch+i] - h_virial_8.data[j*pitch+i]) * double(h_virial_9.data[j*pitch+i] - h_virial_8.data[j*pitch+i]);
 538
 539         // also check that each individual calculation is somewhat close
 540         }
 541     deltaf2 /= double(sysdef->getParticleData()->getN());
 542     deltape2 /= double(sysdef->getParticleData()->getN());
 543     for (unsigned int j = 0; j < 6; j++)
 544         deltav2[j] /= double(sysdef->getParticleData()->getN());
 545     BOOST_CHECK_SMALL(deltaf2, double(tol_small));
 546     BOOST_CHECK_SMALL(deltape2, double(tol_small));
 547     BOOST_CHECK_SMALL(deltav2[0], double(tol_small));
 548     BOOST_CHECK_SMALL(deltav2[1], double(tol_small));
 549     BOOST_CHECK_SMALL(deltav2[2], double(tol_small));
 550     BOOST_CHECK_SMALL(deltav2[3], double(tol_small));
 551     BOOST_CHECK_SMALL(deltav2[4], double(tol_small));
 552     BOOST_CHECK_SMALL(deltav2[5], double(tol_small));
 553     }
 554     }
 555
 556 //! CGCMMForceCompute creator for unit tests
 557 boost::shared_ptr<CGCMMForceCompute> base_class_cgcmm_creator(boost::shared_ptr<SystemDefinition> sysdef, boost::shared_ptr<NeighborList> nlist, Scalar r_cut)
 558     {
 559     return boost::shared_ptr<CGCMMForceCompute>(new CGCMMForceCompute(sysdef, nlist, r_cut));
 560     }
 561
 562 #ifdef ENABLE_CUDA
 563 //! CGCMMForceComputeGPU creator for unit tests
 564 boost::shared_ptr<CGCMMForceCompute> gpu_cgcmm_creator(boost::shared_ptr<SystemDefinition> sysdef, boost::shared_ptr<NeighborList> nlist, Scalar r_cut)
 565     {
 566     nlist->setStorageMode(NeighborList::full);
 567     boost::shared_ptr<CGCMMForceComputeGPU> cgcmm(new CGCMMForceComputeGPU(sysdef, nlist, r_cut));
 568     // the default block size kills valgrind :) reduce it
 569     cgcmm->setBlockSize(64);
 570     return cgcmm;
 571     }
 572 #endif
 573
 574 //! boost test case for particle test on CPU
 575 BOOST_AUTO_TEST_CASE( CGCMMForce_particle124 )
 576     {
 577     cgcmmforce_creator cgcmm_creator_base = bind(base_class_cgcmm_creator, _1, _2, _3);
 578     cgcmm_force_particle124_test(cgcmm_creator_base, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::CPU)));
 579     }
 580
 581 //! boost test case for particle test on CPU
 582 BOOST_AUTO_TEST_CASE( CGCMMForce_particle96 )
 583     {
 584     cgcmmforce_creator cgcmm_creator_base = bind(base_class_cgcmm_creator, _1, _2, _3);
 585     cgcmm_force_particle96_test(cgcmm_creator_base, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::CPU)));
 586     }
 587
 588 //! boost test case for periodic test on CPU
 589 BOOST_AUTO_TEST_CASE( CGCMMForce_periodic )
 590     {
 591     cgcmmforce_creator cgcmm_creator_base = bind(base_class_cgcmm_creator, _1, _2, _3);
 592     cgcmm_force_periodic_test(cgcmm_creator_base, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::CPU)));
 593     }
 594
 595 # ifdef ENABLE_CUDA
 596 //! boost test case for particle test on GPU - threaded
 597 BOOST_AUTO_TEST_CASE( CGCMMForceGPU_particle124 )
 598     {
 599     cgcmmforce_creator cgcmm_creator_gpu = bind(gpu_cgcmm_creator, _1, _2, _3);
 600     cgcmm_force_particle124_test(cgcmm_creator_gpu, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::GPU)));
 601     }
 602
 603 //! boost test case for particle test on GPU - threaded
 604 BOOST_AUTO_TEST_CASE( CGCMMForceGPU_particle96 )
 605     {
 606     cgcmmforce_creator cgcmm_creator_gpu = bind(gpu_cgcmm_creator, _1, _2, _3);
 607     cgcmm_force_particle96_test(cgcmm_creator_gpu, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::GPU)));
 608     }
 609
 610 //! boost test case for periodic test on the GPU
 611 BOOST_AUTO_TEST_CASE( CGCMMForceGPU_periodic )
 612     {
 613     cgcmmforce_creator cgcmm_creator_gpu = bind(gpu_cgcmm_creator, _1, _2, _3);
 614     cgcmm_force_periodic_test(cgcmm_creator_gpu, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::GPU)));
 615     }
 616
 617 //! boost test case for comparing GPU output to base class output
 618 BOOST_AUTO_TEST_CASE( CGCMMForceGPU_compare )
 619     {
 620     cgcmmforce_creator cgcmm_creator_gpu = bind(gpu_cgcmm_creator, _1, _2, _3);
 621     cgcmmforce_creator cgcmm_creator_base = bind(base_class_cgcmm_creator, _1, _2, _3);
 622     cgcmm_force_comparison_test(cgcmm_creator_base, cgcmm_creator_gpu, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::GPU)));
 623     }
 624
 625 #endif
 626
 627 #ifdef WIN32
 628 #pragma warning( pop )
 629 #endif