test/unit/test_gaussian_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 #include <fstream>
  58
  59 #include <boost/bind.hpp>
  60 #include <boost/function.hpp>
  61 #include <boost/shared_ptr.hpp>
  62
  63 #include "AllPairPotentials.h"
  64
  65 #include "NeighborListBinned.h"
  66 #include "Initializers.h"
  67 #include "SnapshotSystemData.h"
  68
  69 #include <math.h>
  70
  71 using namespace std;
  72 using namespace boost;
  73
  74 /*! \file gaussian_force_test.cc
  75     \brief Implements unit tests for PotentialPairGauss and descendants
  76     \ingroup unit_tests
  77 */
  78
  79 //! Name the unit test module
  80 #define BOOST_TEST_MODULE PotentialPairGaussTests
  81 #include "boost_utf_configure.h"
  82
  83 //! Typedef'd PotentialPairGauss factory
  84 typedef boost::function<boost::shared_ptr<PotentialPairGauss> (boost::shared_ptr<SystemDefinition> sysdef,
  85                                                         boost::shared_ptr<NeighborList> nlist)> gaussforce_creator;
  86
  87 //! Test the ability of the gauss force compute to actually calucate forces
  88 void gauss_force_particle_test(gaussforce_creator gauss_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     pdata_3->setFlags(~PDataFlags(0));
 101
 102     pdata_3->setPosition(0,make_scalar3(0.0,0.0,0.0));
 103     pdata_3->setPosition(1,make_scalar3(1.0,0.0,0.0));
 104     pdata_3->setPosition(2,make_scalar3(2.0,0.0,0.0));
 105
 106     boost::shared_ptr<NeighborList> nlist_3(new NeighborList(sysdef_3, Scalar(1.3), Scalar(3.0)));
 107     boost::shared_ptr<PotentialPairGauss> fc_3 = gauss_creator(sysdef_3, nlist_3);
 108     fc_3->setRcut(0, 0, Scalar(1.3));
 109
 110     // first test: choose a basic sigma
 111     Scalar epsilon = Scalar(1.15);
 112     Scalar sigma = Scalar(0.5);
 113     fc_3->setParams(0,0,make_scalar2(epsilon,sigma));
 114
 115     // compute the forces
 116     fc_3->compute(0);
 117
 118     {
 119     GPUArray<Scalar4>& force_array_1 =  fc_3->getForceArray();
 120     GPUArray<Scalar>& virial_array_1 =  fc_3->getVirialArray();
 121     unsigned int pitch = virial_array_1.getPitch();
 122     ArrayHandle<Scalar4> h_force_1(force_array_1,access_location::host,access_mode::read);
 123     ArrayHandle<Scalar> h_virial_1(virial_array_1,access_location::host,access_mode::read);
 124     MY_BOOST_CHECK_CLOSE(h_force_1.data[0].x, -0.622542302888418, tol);
 125     MY_BOOST_CHECK_SMALL(h_force_1.data[0].y, tol_small);
 126     MY_BOOST_CHECK_SMALL(h_force_1.data[0].z, tol_small);
 127     MY_BOOST_CHECK_CLOSE(h_force_1.data[0].w, 0.155635575722105/2.0, tol);
 128     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_1.data[0*pitch+0]
 129                                        +h_virial_1.data[3*pitch+0]
 130                                        +h_virial_1.data[5*pitch+0]), 0.103757050481403, tol);
 131
 132     MY_BOOST_CHECK_SMALL(h_force_1.data[1].x, tol_small);
 133     MY_BOOST_CHECK_SMALL(h_force_1.data[1].y, tol_small);
 134     MY_BOOST_CHECK_SMALL(h_force_1.data[1].z, tol_small);
 135     MY_BOOST_CHECK_CLOSE(h_force_1.data[1].w, 0.155635575722105, tol);
 136     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_1.data[0*pitch+1]
 137                                        +h_virial_1.data[3*pitch+1]
 138                                        +h_virial_1.data[5*pitch+1]), 0.103757050481403*2, tol);
 139
 140     MY_BOOST_CHECK_CLOSE(h_force_1.data[2].x, 0.622542302888418, tol);
 141     MY_BOOST_CHECK_SMALL(h_force_1.data[2].y, tol_small);
 142     MY_BOOST_CHECK_SMALL(h_force_1.data[2].z, tol_small);
 143     MY_BOOST_CHECK_CLOSE(h_force_1.data[2].w, 0.155635575722105/2.0, tol);
 144     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_1.data[0*pitch+2]
 145                                        +h_virial_1.data[3*pitch+2]
 146                                        +h_virial_1.data[5*pitch+2]), 0.103757050481403, tol);
 147     }
 148
 149     // swap the order of particles 0 ans 2 in memory to check that the force compute handles this properly
 150     {
 151     ArrayHandle<Scalar4> h_pos(pdata_3->getPositions(), access_location::host, access_mode::readwrite);
 152     ArrayHandle<unsigned int> h_tag(pdata_3->getTags(), access_location::host, access_mode::readwrite);
 153     ArrayHandle<unsigned int> h_rtag(pdata_3->getRTags(), access_location::host, access_mode::readwrite);
 154
 155     h_pos.data[2].x = h_pos.data[2].y = h_pos.data[2].z = 0.0;
 156     h_pos.data[0].x = Scalar(2.0); h_pos.data[0].y = h_pos.data[0].z = 0.0;
 157
 158     h_tag.data[0] = 2;
 159     h_tag.data[2] = 0;
 160     h_rtag.data[0] = 2;
 161     h_rtag.data[2] = 0;
 162     }
 163
 164     // notify the particle data that we changed the order
 165     pdata_3->notifyParticleSort();
 166
 167     // recompute the forces at the same timestep, they should be updated
 168     fc_3->compute(1);
 169
 170     {
 171     GPUArray<Scalar4>& force_array_2 =  fc_3->getForceArray();
 172     GPUArray<Scalar>& virial_array_2 =  fc_3->getVirialArray();
 173     ArrayHandle<Scalar4> h_force_2(force_array_2,access_location::host,access_mode::read);
 174     ArrayHandle<Scalar> h_virial_2(virial_array_2,access_location::host,access_mode::read);
 175     MY_BOOST_CHECK_CLOSE(h_force_2.data[0].x, 0.622542302888418, tol);
 176     MY_BOOST_CHECK_CLOSE(h_force_2.data[2].x, -0.622542302888418, tol);
 177     }
 178     }
 179
 180 //! Tests the ability of a PotentialPairGauss to handle periodic boundary conditions
 181 void gauss_force_periodic_test(gaussforce_creator gauss_creator, boost::shared_ptr<ExecutionConfiguration> exec_conf)
 182     {
 183     ////////////////////////////////////////////////////////////////////
 184     // now, lets do a more thorough test and include boundary conditions
 185     // there are way too many permutations to test here, so I will simply
 186     // test +x, -x, +y, -y, +z, and -z independantly
 187     // build a 6 particle system with particles across each boundary
 188     // also test the ability of the force compute to use different particle types
 189     boost::shared_ptr<SystemDefinition> sysdef_6(new SystemDefinition(6, BoxDim(20.0, 40.0, 60.0), 3, 0, 0, 0, 0, exec_conf));
 190     boost::shared_ptr<ParticleData> pdata_6 = sysdef_6->getParticleData();
 191     pdata_6->setFlags(~PDataFlags(0));
 192
 193     pdata_6->setPosition(0, make_scalar3(-9.6,0.0,0.0));
 194     pdata_6->setPosition(1, make_scalar3(9.6, 0.0,0.0));
 195     pdata_6->setPosition(2, make_scalar3(0.0,-19.6,0.0));
 196     pdata_6->setPosition(3, make_scalar3(0.0,19.6,0.0));
 197     pdata_6->setPosition(4, make_scalar3(0.0,0.0,-29.6));
 198     pdata_6->setPosition(5, make_scalar3(0.0,0.0,29.6));
 199
 200     pdata_6->setType(0,0);
 201     pdata_6->setType(1,1);
 202     pdata_6->setType(2,2);
 203     pdata_6->setType(3,0);
 204     pdata_6->setType(4,2);
 205     pdata_6->setType(5,1);
 206
 207     boost::shared_ptr<NeighborList> nlist_6(new NeighborList(sysdef_6, Scalar(1.3), Scalar(3.0)));
 208     boost::shared_ptr<PotentialPairGauss> fc_6 = gauss_creator(sysdef_6, nlist_6);
 209     fc_6->setRcut(0, 0, Scalar(1.3));
 210     fc_6->setRcut(0, 1, Scalar(1.3));
 211     fc_6->setRcut(0, 2, Scalar(1.3));
 212     fc_6->setRcut(1, 1, Scalar(1.3));
 213     fc_6->setRcut(1, 2, Scalar(1.3));
 214     fc_6->setRcut(2, 2, Scalar(1.3));
 215
 216     // choose a small sigma so that all interactions are attractive
 217     Scalar epsilon = Scalar(1.0);
 218     Scalar sigma = Scalar(0.5);
 219
 220     // make life easy: just change epsilon for the different pairs
 221     fc_6->setParams(0,0,make_scalar2(epsilon,sigma));
 222     fc_6->setParams(0,1,make_scalar2(Scalar(2.0)*epsilon,sigma));
 223     fc_6->setParams(0,2,make_scalar2(Scalar(3.0)*epsilon,sigma));
 224     fc_6->setParams(1,1,make_scalar2(Scalar(4.0)*epsilon,sigma));
 225     fc_6->setParams(1,2,make_scalar2(Scalar(5.0)*epsilon,sigma));
 226     fc_6->setParams(2,2,make_scalar2(Scalar(6.0)*epsilon,sigma));
 227
 228     fc_6->compute(0);
 229
 230     {
 231     GPUArray<Scalar4>& force_array_3 =  fc_6->getForceArray();
 232     GPUArray<Scalar>& virial_array_3 =  fc_6->getVirialArray();
 233     unsigned int pitch = virial_array_3.getPitch();
 234     ArrayHandle<Scalar4> h_force_3(force_array_3,access_location::host,access_mode::read);
 235     ArrayHandle<Scalar> h_virial_3(virial_array_3,access_location::host,access_mode::read);
 236     // particle 0 should be pushed right
 237     MY_BOOST_CHECK_CLOSE(h_force_3.data[0].x, 2.224298403625553*0.8, tol);
 238     MY_BOOST_CHECK_SMALL(h_force_3.data[0].y, tol_small);
 239     MY_BOOST_CHECK_SMALL(h_force_3.data[0].z, tol_small);
 240     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_3.data[0*pitch+0]
 241                                        +h_virial_3.data[3*pitch+0]
 242                                        +h_virial_3.data[5*pitch+0]), 0.296573120483407*0.8, tol);
 243
 244     // particle 1 should be pushed left
 245     MY_BOOST_CHECK_CLOSE(h_force_3.data[1].x, -2.224298403625553*0.8, tol);
 246     MY_BOOST_CHECK_SMALL(h_force_3.data[1].y, tol_small);
 247     MY_BOOST_CHECK_SMALL(h_force_3.data[1].z, tol_small);
 248     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_3.data[0*pitch+1]
 249                                        +h_virial_3.data[3*pitch+1]
 250                                        +h_virial_3.data[5*pitch+1]), 0.296573120483407*0.8, tol);
 251
 252     // particle 2 should be pushed up
 253     MY_BOOST_CHECK_CLOSE(h_force_3.data[2].y, 3.336447605438329*0.8, tol);
 254     MY_BOOST_CHECK_SMALL(h_force_3.data[2].x, tol_small);
 255     MY_BOOST_CHECK_SMALL(h_force_3.data[2].z, tol_small);
 256     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_3.data[0*pitch+2]
 257                                        +h_virial_3.data[3*pitch+2]
 258                                        +h_virial_3.data[5*pitch+2]), 0.444859680725111*0.8, tol);
 259
 260     // particle 3 should be pushed down
 261     MY_BOOST_CHECK_CLOSE(h_force_3.data[3].y, -3.336447605438329*0.8, tol);
 262     MY_BOOST_CHECK_SMALL(h_force_3.data[3].x, tol_small);
 263     MY_BOOST_CHECK_SMALL(h_force_3.data[3].z, tol_small);
 264     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_3.data[0*pitch+3]
 265                                        +h_virial_3.data[3*pitch+3]
 266                                        +h_virial_3.data[5*pitch+3]), 0.444859680725111*0.8, tol);
 267
 268     // particle 4 should be pushed forward
 269     MY_BOOST_CHECK_CLOSE(h_force_3.data[4].z, 5.560746009063882*0.8, tol);
 270     MY_BOOST_CHECK_SMALL(h_force_3.data[4].x, tol_small);
 271     MY_BOOST_CHECK_SMALL(h_force_3.data[4].y, tol_small);
 272     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_3.data[0*pitch+4]
 273                                        +h_virial_3.data[3*pitch+4]
 274                                        +h_virial_3.data[5*pitch+4]),  0.741432801208518*0.8, tol);
 275
 276     // particle 3 should be pushed back
 277     MY_BOOST_CHECK_CLOSE(h_force_3.data[5].z, -5.560746009063882*0.8, tol);
 278     MY_BOOST_CHECK_SMALL(h_force_3.data[5].x, tol_small);
 279     MY_BOOST_CHECK_SMALL(h_force_3.data[5].y, tol_small);
 280     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_3.data[0*pitch+5]
 281                                        +h_virial_3.data[3*pitch+5]
 282                                        +h_virial_3.data[5*pitch+5]),  0.741432801208518*0.8, tol);
 283     }
 284     }
 285
 286 //! Unit test a comparison between 2 LJForceComputes on a "real" system
 287 void gauss_force_comparison_test(gaussforce_creator gauss_creator1,
 288                                  gaussforce_creator gauss_creator2,
 289                                  boost::shared_ptr<ExecutionConfiguration> exec_conf)
 290     {
 291     const unsigned int N = 5000;
 292
 293     // create a random particle system to sum forces on
 294     RandomInitializer rand_init(N, Scalar(0.2), Scalar(0.9), "A");
 295     boost::shared_ptr<SnapshotSystemData> snap;
 296     snap = rand_init.getSnapshot();
 297     boost::shared_ptr<SystemDefinition> sysdef(new SystemDefinition(snap, exec_conf));
 298     boost::shared_ptr<ParticleData> pdata = sysdef->getParticleData();
 299     pdata->setFlags(~PDataFlags(0));
 300     boost::shared_ptr<NeighborListBinned> nlist(new NeighborListBinned(sysdef, Scalar(3.0), Scalar(0.8)));
 301
 302     boost::shared_ptr<PotentialPairGauss> fc1 = gauss_creator1(sysdef, nlist);
 303     boost::shared_ptr<PotentialPairGauss> fc2 = gauss_creator2(sysdef, nlist);
 304     fc1->setRcut(0, 0, Scalar(3.0));
 305     fc2->setRcut(0, 0, Scalar(3.0));
 306
 307     // setup some values for epsilon and sigma
 308     Scalar epsilon = Scalar(1.0);
 309     Scalar sigma = Scalar(1.2);
 310
 311     // specify the force parameters
 312     fc1->setParams(0,0,make_scalar2(epsilon,sigma));
 313     fc2->setParams(0,0,make_scalar2(epsilon,sigma));
 314
 315     // compute the forces
 316     fc1->compute(0);
 317     fc2->compute(0);
 318
 319     {
 320     // verify that the forces are identical (within roundoff errors)
 321     GPUArray<Scalar4>& force_array_4 =  fc1->getForceArray();
 322     GPUArray<Scalar>& virial_array_4 =  fc1->getVirialArray();
 323     unsigned int pitch = virial_array_4.getPitch();
 324     ArrayHandle<Scalar4> h_force_4(force_array_4,access_location::host,access_mode::read);
 325     ArrayHandle<Scalar> h_virial_4(virial_array_4,access_location::host,access_mode::read);
 326     GPUArray<Scalar4>& force_array_5 =  fc2->getForceArray();
 327     GPUArray<Scalar>& virial_array_5 =  fc2->getVirialArray();
 328     ArrayHandle<Scalar4> h_force_5(force_array_5,access_location::host,access_mode::read);
 329     ArrayHandle<Scalar> h_virial_5(virial_array_5,access_location::host,access_mode::read);
 330
 331     // compare average deviation between the two computes
 332     double deltaf2 = 0.0;
 333     double deltape2 = 0.0;
 334     double deltav2[6];
 335     for (unsigned int i = 0; i < 6;  i++)
 336         deltav2[i] = 0.0;
 337
 338     for (unsigned int i = 0; i < N; i++)
 339         {
 340         deltaf2 += double(h_force_5.data[i].x - h_force_4.data[i].x) * double(h_force_5.data[i].x - h_force_4.data[i].x);
 341         deltaf2 += double(h_force_5.data[i].y - h_force_4.data[i].y) * double(h_force_5.data[i].y - h_force_4.data[i].y);
 342         deltaf2 += double(h_force_5.data[i].z - h_force_4.data[i].z) * double(h_force_5.data[i].z - h_force_4.data[i].z);
 343         deltape2 += double(h_force_5.data[i].w - h_force_4.data[i].w) * double(h_force_5.data[i].w - h_force_4.data[i].w);
 344         for (unsigned int j = 0; j < 6; j++)
 345             deltav2[j] += double(h_virial_5.data[j*pitch+i] - h_virial_4.data[j*pitch+i]) * double(h_virial_5.data[j*pitch+i] - h_virial_4.data[j*pitch+i]);
 346
 347         // also check that each individual calculation is somewhat close
 348         }
 349     deltaf2 /= double(pdata->getN());
 350     deltape2 /= double(pdata->getN());
 351     for (unsigned int i = 0; i < 6; i++)
 352         deltav2[i] /= double(pdata->getN());
 353     BOOST_CHECK_SMALL(deltaf2, double(tol_small));
 354     BOOST_CHECK_SMALL(deltape2, double(tol_small));
 355     BOOST_CHECK_SMALL(deltav2[0], double(tol_small));
 356     BOOST_CHECK_SMALL(deltav2[1], double(tol_small));
 357     BOOST_CHECK_SMALL(deltav2[2], double(tol_small));
 358     BOOST_CHECK_SMALL(deltav2[3], double(tol_small));
 359     BOOST_CHECK_SMALL(deltav2[4], double(tol_small));
 360     BOOST_CHECK_SMALL(deltav2[5], double(tol_small));
 361     }
 362     }
 363
 364 //! Test the ability of the gauss force compute to compute forces with different shift modes
 365 void gauss_force_shift_test(gaussforce_creator gauss_creator, boost::shared_ptr<ExecutionConfiguration> exec_conf)
 366     {
 367     // this 2-particle test is just to get a plot of the potential and force vs r cut
 368     boost::shared_ptr<SystemDefinition> sysdef_2(new SystemDefinition(2, BoxDim(1000.0), 1, 0, 0, 0, 0, exec_conf));
 369     boost::shared_ptr<ParticleData> pdata_2 = sysdef_2->getParticleData();
 370     pdata_2->setFlags(~PDataFlags(0));
 371
 372     pdata_2->setPosition(0,make_scalar3(0.0,0.0,0.0));
 373     pdata_2->setPosition(1,make_scalar3(2.8,0.0,0.0));
 374     boost::shared_ptr<NeighborList> nlist_2(new NeighborList(sysdef_2, Scalar(3.0), Scalar(0.8)));
 375     boost::shared_ptr<PotentialPairGauss> fc_no_shift = gauss_creator(sysdef_2, nlist_2);
 376     fc_no_shift->setShiftMode(PotentialPairGauss::no_shift);
 377     fc_no_shift->setRcut(0, 0, Scalar(3.0));
 378     boost::shared_ptr<PotentialPairGauss> fc_shift = gauss_creator(sysdef_2, nlist_2);
 379     fc_shift->setShiftMode(PotentialPairGauss::shift);
 380     fc_shift->setRcut(0, 0, Scalar(3.0));
 381
 382     nlist_2->setStorageMode(NeighborList::full);
 383
 384     // setup a standard epsilon and sigma
 385     Scalar epsilon = Scalar(1.0);
 386     Scalar sigma = Scalar(1.0);
 387     fc_no_shift->setParams(0,0,make_scalar2(epsilon,sigma));
 388     fc_shift->setParams(0,0,make_scalar2(epsilon,sigma));
 389
 390     fc_no_shift->compute(0);
 391     fc_shift->compute(0);
 392
 393     {
 394     GPUArray<Scalar4>& force_array_6 =  fc_no_shift->getForceArray();
 395     GPUArray<Scalar>& virial_array_6 =  fc_no_shift->getVirialArray();
 396     ArrayHandle<Scalar4> h_force_6(force_array_6,access_location::host,access_mode::read);
 397     ArrayHandle<Scalar> h_virial_6(virial_array_6,access_location::host,access_mode::read);
 398
 399     MY_BOOST_CHECK_CLOSE(h_force_6.data[0].x, -0.055555065284237, tol);
 400     MY_BOOST_CHECK_CLOSE(h_force_6.data[0].w, 0.019841094744370/2.0, tol);
 401     MY_BOOST_CHECK_CLOSE(h_force_6.data[1].x, 0.055555065284237, tol);
 402     MY_BOOST_CHECK_CLOSE(h_force_6.data[1].w, 0.019841094744370/2.0, tol);
 403
 404     GPUArray<Scalar4>& force_array_7 =  fc_shift->getForceArray();
 405     GPUArray<Scalar>& virial_array_7 =  fc_shift->getVirialArray();
 406     ArrayHandle<Scalar4> h_force_7(force_array_7,access_location::host,access_mode::read);
 407     ArrayHandle<Scalar> h_virial_7(virial_array_7,access_location::host,access_mode::read);
 408
 409     // shifted just has pe shifted by a given amount
 410     MY_BOOST_CHECK_CLOSE(h_force_7.data[0].x, -0.055555065284237, tol);
 411     MY_BOOST_CHECK_CLOSE(h_force_7.data[0].w, 0.008732098206128/2.0, tol);
 412     MY_BOOST_CHECK_CLOSE(h_force_7.data[1].x, 0.055555065284237, tol);
 413     MY_BOOST_CHECK_CLOSE(h_force_7.data[1].w, 0.008732098206128/2.0, tol);
 414     }
 415
 416     // check once again to verify that nothing fish happens past r_cut
 417     pdata_2->setPosition(0,make_scalar3(0.0,0.0,0.0));
 418     pdata_2->setPosition(1,make_scalar3(3.1,0.0,0.0));
 419
 420     fc_no_shift->compute(2);
 421     fc_shift->compute(2);
 422
 423     {
 424     GPUArray<Scalar4>& force_array_8 =  fc_no_shift->getForceArray();
 425     GPUArray<Scalar>& virial_array_8 =  fc_no_shift->getVirialArray();
 426     ArrayHandle<Scalar4> h_force_8(force_array_8,access_location::host,access_mode::read);
 427     ArrayHandle<Scalar> h_virial_8(virial_array_8,access_location::host,access_mode::read);
 428     GPUArray<Scalar4>& force_array_9 =  fc_shift->getForceArray();
 429     GPUArray<Scalar>& virial_array_9 =  fc_shift->getVirialArray();
 430     ArrayHandle<Scalar4> h_force_9(force_array_9,access_location::host,access_mode::read);
 431     ArrayHandle<Scalar> h_virial_9(virial_array_9,access_location::host,access_mode::read);
 432
 433     MY_BOOST_CHECK_SMALL(h_force_9.data[0].x, tol_small);
 434     MY_BOOST_CHECK_SMALL(h_force_9.data[0].w, tol_small);
 435     MY_BOOST_CHECK_SMALL(h_force_9.data[1].x, tol_small);
 436     MY_BOOST_CHECK_SMALL(h_force_9.data[1].w, tol_small);
 437
 438     // shifted just has pe shifted by a given amount
 439     MY_BOOST_CHECK_SMALL(h_force_9.data[0].x, tol_small);
 440     MY_BOOST_CHECK_SMALL(h_force_9.data[0].w, tol_small);
 441     MY_BOOST_CHECK_SMALL(h_force_9.data[1].x, tol_small);
 442     MY_BOOST_CHECK_SMALL(h_force_9.data[1].w, tol_small);
 443     }
 444     }
 445
 446 //! LJForceCompute creator for unit tests
 447 boost::shared_ptr<PotentialPairGauss> base_class_gauss_creator(boost::shared_ptr<SystemDefinition> sysdef,
 448                                                         boost::shared_ptr<NeighborList> nlist)
 449     {
 450     return boost::shared_ptr<PotentialPairGauss>(new PotentialPairGauss(sysdef, nlist));
 451     }
 452
 453 #ifdef ENABLE_CUDA
 454 //! PotentialPairGaussGPU creator for unit tests
 455 boost::shared_ptr<PotentialPairGaussGPU> gpu_gauss_creator(boost::shared_ptr<SystemDefinition> sysdef,
 456                                                     boost::shared_ptr<NeighborList> nlist)
 457     {
 458     nlist->setStorageMode(NeighborList::full);
 459     boost::shared_ptr<PotentialPairGaussGPU> gauss(new PotentialPairGaussGPU(sysdef, nlist));
 460     return gauss;
 461     }
 462 #endif
 463
 464 //! boost test case for particle test on CPU
 465 BOOST_AUTO_TEST_CASE( GaussForce_particle )
 466     {
 467     gaussforce_creator gauss_creator_base = bind(base_class_gauss_creator, _1, _2);
 468     gauss_force_particle_test(gauss_creator_base, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::CPU)));
 469     }
 470
 471 //! boost test case for periodic test on CPU
 472 BOOST_AUTO_TEST_CASE( GaussForce_periodic )
 473     {
 474     gaussforce_creator gauss_creator_base = bind(base_class_gauss_creator, _1, _2);
 475     gauss_force_periodic_test(gauss_creator_base, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::CPU)));
 476     }
 477
 478 //! boost test case for particle test on CPU
 479 BOOST_AUTO_TEST_CASE( GaussForce_shift )
 480     {
 481     gaussforce_creator gauss_creator_base = bind(base_class_gauss_creator, _1, _2);
 482     gauss_force_shift_test(gauss_creator_base, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::CPU)));
 483     }
 484
 485 # ifdef ENABLE_CUDA
 486 //! boost test case for particle test on GPU
 487 BOOST_AUTO_TEST_CASE( GaussForceGPU_particle )
 488     {
 489     gaussforce_creator gauss_creator_gpu = bind(gpu_gauss_creator, _1, _2);
 490     gauss_force_particle_test(gauss_creator_gpu, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::GPU)));
 491     }
 492
 493 //! boost test case for periodic test on the GPU
 494 BOOST_AUTO_TEST_CASE( GaussForceGPU_periodic )
 495     {
 496     gaussforce_creator gauss_creator_gpu = bind(gpu_gauss_creator, _1, _2);
 497     gauss_force_periodic_test(gauss_creator_gpu, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::GPU)));
 498     }
 499
 500 //! boost test case for shift test on GPU
 501 BOOST_AUTO_TEST_CASE( GaussForceGPU_shift )
 502     {
 503     gaussforce_creator gauss_creator_gpu = bind(gpu_gauss_creator, _1, _2);
 504     gauss_force_shift_test(gauss_creator_gpu, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::GPU)));
 505     }
 506
 507 //! boost test case for comparing GPU output to base class output
 508 BOOST_AUTO_TEST_CASE( GaussForceGPU_compare )
 509     {
 510     gaussforce_creator gauss_creator_gpu = bind(gpu_gauss_creator, _1, _2);
 511     gaussforce_creator gauss_creator_base = bind(base_class_gauss_creator, _1, _2);
 512     gauss_force_comparison_test(gauss_creator_base, gauss_creator_gpu, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::GPU)));
 513     }
 514
 515 #endif
 516
 517 #ifdef WIN32
 518 #pragma warning( pop )
 519 #endif