test/unit/test_harmonic_bond_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
  61 #include "AllBondPotentials.h"
  62 #include "ConstForceCompute.h"
  63 #include "SnapshotSystemData.h"
  64
  65 #include "Initializers.h"
  66
  67 using namespace std;
  68 using namespace boost;
  69
  70 /*! \file harmonic_bond_force_test.cc
  71     \brief Implements unit tests for PotentialBondHarmonic and
  72            PotentialBondHarmonicGPU
  73     \ingroup unit_tests
  74 */
  75
  76 //! Name the boost unit test module
  77 #define BOOST_TEST_MODULE BondForceTests
  78 #include "boost_utf_configure.h"
  79
  80 //! Typedef to make using the boost::function factory easier
  81 typedef boost::function<boost::shared_ptr<PotentialBondHarmonic>  (boost::shared_ptr<SystemDefinition> sysdef)> bondforce_creator;
  82
  83 //! Perform some simple functionality tests of any BondForceCompute
  84 void bond_force_basic_tests(bondforce_creator bf_creator, boost::shared_ptr<ExecutionConfiguration> exec_conf)
  85     {
  86     /////////////////////////////////////////////////////////
  87     // start with the simplest possible test: 2 particles in a huge box with only one bond type
  88     boost::shared_ptr<SystemDefinition> sysdef_2(new SystemDefinition(2, BoxDim(1000.0), 1, 1, 0, 0, 0, exec_conf));
  89     boost::shared_ptr<ParticleData> pdata_2 = sysdef_2->getParticleData();
  90
  91     pdata_2->setPosition(0,make_scalar3(0.0,0.0,0.0));
  92     pdata_2->setPosition(1,make_scalar3(0.9,0.0,0.0));
  93     pdata_2->setFlags(~PDataFlags(0));
  94
  95     // create the bond force compute to check
  96     boost::shared_ptr<PotentialBondHarmonic> fc_2 = bf_creator(sysdef_2);
  97     fc_2->setParams(0, make_scalar2(1.5, 0.75));
  98
  99     // compute the force and check the results
 100     fc_2->compute(0);
 101     GPUArray<Scalar4>& force_array_1 =  fc_2->getForceArray();
 102     GPUArray<Scalar>& virial_array_1 =  fc_2->getVirialArray();
 103
 104     {
 105     unsigned int pitch = virial_array_1.getPitch();
 106     ArrayHandle<Scalar4> h_force_1(force_array_1,access_location::host,access_mode::read);
 107     ArrayHandle<Scalar> h_virial_1(virial_array_1,access_location::host,access_mode::read);
 108     // check that the force is correct, it should be 0 since we haven't created any bonds yet
 109     MY_BOOST_CHECK_SMALL(h_force_1.data[0].x, tol_small);
 110     MY_BOOST_CHECK_SMALL(h_force_1.data[0].y, tol_small);
 111     MY_BOOST_CHECK_SMALL(h_force_1.data[0].z, tol_small);
 112     MY_BOOST_CHECK_SMALL(h_force_1.data[0].w, tol_small);
 113     MY_BOOST_CHECK_SMALL(h_virial_1.data[0*pitch+0], tol_small);
 114     MY_BOOST_CHECK_SMALL(h_virial_1.data[1*pitch+0], tol_small);
 115     MY_BOOST_CHECK_SMALL(h_virial_1.data[2*pitch+0], tol_small);
 116     MY_BOOST_CHECK_SMALL(h_virial_1.data[3*pitch+0], tol_small);
 117     MY_BOOST_CHECK_SMALL(h_virial_1.data[4*pitch+0], tol_small);
 118     MY_BOOST_CHECK_SMALL(h_virial_1.data[5*pitch+0], tol_small);
 119     }
 120
 121     // add a bond and check again
 122     sysdef_2->getBondData()->addBondedGroup(Bond(0, 0,1));
 123     fc_2->compute(1);
 124
 125     {
 126     // this time there should be a force
 127     GPUArray<Scalar4>& force_array_2 =  fc_2->getForceArray();
 128     GPUArray<Scalar>& virial_array_2 =  fc_2->getVirialArray();
 129     unsigned int pitch = virial_array_2.getPitch();
 130     ArrayHandle<Scalar4> h_force_2(force_array_2,access_location::host,access_mode::read);
 131     ArrayHandle<Scalar> h_virial_2(virial_array_2,access_location::host,access_mode::read);
 132     MY_BOOST_CHECK_CLOSE(h_force_2.data[0].x, 0.225, tol);
 133     MY_BOOST_CHECK_SMALL(h_force_2.data[0].y, tol_small);
 134     MY_BOOST_CHECK_SMALL(h_force_2.data[0].z, tol_small);
 135     MY_BOOST_CHECK_CLOSE(h_force_2.data[0].w, 0.0084375, tol);
 136     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_2.data[0*pitch+0]
 137                                        +h_virial_2.data[3*pitch+0]
 138                                        +h_virial_2.data[5*pitch+0]), -0.03375, tol);
 139
 140     // check that the two forces are negatives of each other
 141     MY_BOOST_CHECK_CLOSE(h_force_2.data[0].x, -h_force_2.data[1].x, tol);
 142     MY_BOOST_CHECK_CLOSE(h_force_2.data[0].y, -h_force_2.data[1].y, tol);
 143     MY_BOOST_CHECK_CLOSE(h_force_2.data[0].z, -h_force_2.data[1].z, tol);
 144     MY_BOOST_CHECK_CLOSE(h_force_2.data[0].w, h_force_2.data[1].w, tol);
 145     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_2.data[0*pitch+1]
 146                                        +h_virial_2.data[3*pitch+1]
 147                                        +h_virial_2.data[5*pitch+1]), -0.03375, tol);
 148     }
 149
 150     // rearrange the two particles in memory and see if they are properly updated
 151     {
 152     ArrayHandle<Scalar4> h_pos(pdata_2->getPositions(), access_location::host, access_mode::readwrite);
 153     ArrayHandle<unsigned int> h_tag(pdata_2->getTags(), access_location::host, access_mode::readwrite);
 154     ArrayHandle<unsigned int> h_rtag(pdata_2->getRTags(), access_location::host, access_mode::readwrite);
 155
 156     h_pos.data[0].x = Scalar(0.9);
 157     h_pos.data[1].x = Scalar(0.0);
 158     h_tag.data[0] = 1;
 159     h_tag.data[1] = 0;
 160     h_rtag.data[0] = 1;
 161     h_rtag.data[1] = 0;
 162     }
 163
 164     // notify that we made the sort
 165     pdata_2->notifyParticleSort();
 166     // recompute at the same timestep, the forces should still be updated
 167     fc_2->compute(1);
 168
 169     {
 170     // this time there should be a force
 171     GPUArray<Scalar4>& force_array_3 =  fc_2->getForceArray();
 172     GPUArray<Scalar>& virial_array_3 =  fc_2->getVirialArray();
 173     ArrayHandle<Scalar4> h_force_3(force_array_3,access_location::host,access_mode::read);
 174     ArrayHandle<Scalar> h_virial_3(virial_array_3,access_location::host,access_mode::read);
 175     MY_BOOST_CHECK_CLOSE(h_force_3.data[0].x, -0.225, tol);
 176     MY_BOOST_CHECK_CLOSE(h_force_3.data[1].x, 0.225, tol);
 177     }
 178
 179     // check r=r_0 behavior
 180     pdata_2->setPosition(0,make_scalar3(0.0,0.0,0.0));
 181     pdata_2->setPosition(1,make_scalar3(0.75,0.0,0.0));
 182
 183     fc_2->compute(2);
 184
 185     {
 186     // the force should be zero
 187     GPUArray<Scalar4>& force_array_4 =  fc_2->getForceArray();
 188     GPUArray<Scalar>& virial_array_4 =  fc_2->getVirialArray();
 189     ArrayHandle<Scalar4> h_force_4(force_array_4,access_location::host,access_mode::read);
 190     ArrayHandle<Scalar> h_virial_4(virial_array_4,access_location::host,access_mode::read);
 191     MY_BOOST_CHECK_SMALL(h_force_4.data[0].x, tol_small);
 192     MY_BOOST_CHECK_SMALL(h_force_4.data[1].x, tol_small);
 193     }
 194
 195     ////////////////////////////////////////////////////////////////////
 196     // now, lets do a more thorough test and include boundary conditions
 197     // there are way too many permutations to test here, so I will simply
 198     // test +x, -x, +y, -y, +z, and -z independantly
 199     // build a 6 particle system with particles across each boundary
 200     // also test more than one type of bond
 201     boost::shared_ptr<SystemDefinition> sysdef_6(new SystemDefinition(6, BoxDim(20.0, 40.0, 60.0), 1, 3, 0, 0, 0, exec_conf));
 202     boost::shared_ptr<ParticleData> pdata_6 = sysdef_6->getParticleData();
 203     pdata_6->setFlags(~PDataFlags(0));
 204
 205     pdata_6->setPosition(0, make_scalar3(-9.6,0.0,0.0));
 206     pdata_6->setPosition(1, make_scalar3(9.6, 0.0,0.0));
 207     pdata_6->setPosition(2, make_scalar3(0.0,-19.6,0.0));
 208     pdata_6->setPosition(3, make_scalar3(0.0,19.6,0.0));
 209     pdata_6->setPosition(4, make_scalar3(0.0,0.0,-29.6));
 210     pdata_6->setPosition(5, make_scalar3(0.0,0.0,29.6));
 211
 212     boost::shared_ptr<PotentialBondHarmonic> fc_6 = bf_creator(sysdef_6);
 213     fc_6->setParams(0, make_scalar2( 1.5, 0.75));
 214     fc_6->setParams(1, make_scalar2(2.0*1.5, 0.75));
 215     fc_6->setParams(2, make_scalar2(1.5, 0.5));
 216
 217     sysdef_6->getBondData()->addBondedGroup(Bond(0, 0,1));
 218     sysdef_6->getBondData()->addBondedGroup(Bond(1, 2,3));
 219     sysdef_6->getBondData()->addBondedGroup(Bond(2, 4,5));
 220
 221     fc_6->compute(0);
 222
 223     {
 224     // check that the forces are correctly computed
 225     GPUArray<Scalar4>& force_array_5 =  fc_6->getForceArray();
 226     GPUArray<Scalar>& virial_array_5 =  fc_6->getVirialArray();
 227     unsigned int pitch = virial_array_5.getPitch();
 228     ArrayHandle<Scalar4> h_force_5(force_array_5,access_location::host,access_mode::read);
 229     ArrayHandle<Scalar> h_virial_5(virial_array_5,access_location::host,access_mode::read);
 230     MY_BOOST_CHECK_CLOSE(h_force_5.data[0].x, -0.075, tol);
 231     MY_BOOST_CHECK_SMALL(h_force_5.data[0].y, tol_small);
 232     MY_BOOST_CHECK_SMALL(h_force_5.data[0].z, tol_small);
 233     MY_BOOST_CHECK_CLOSE(h_force_5.data[0].w, 9.375e-4, tol);
 234     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_5.data[0*pitch+0]
 235                                        +h_virial_5.data[3*pitch+0]
 236                                        +h_virial_5.data[5*pitch+0]), -0.01, tol);
 237
 238     MY_BOOST_CHECK_CLOSE(h_force_5.data[1].x, 0.075, tol);
 239     MY_BOOST_CHECK_SMALL(h_force_5.data[1].y, tol_small);
 240     MY_BOOST_CHECK_SMALL(h_force_5.data[1].z, tol_small);
 241     MY_BOOST_CHECK_CLOSE(h_force_5.data[1].w, 9.375e-4, tol);
 242     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_5.data[0*pitch+1]
 243                                        +h_virial_5.data[3*pitch+1]
 244                                        +h_virial_5.data[5*pitch+1]), -0.01, tol);
 245
 246     MY_BOOST_CHECK_SMALL(h_force_5.data[2].x, tol_small);
 247     MY_BOOST_CHECK_CLOSE(h_force_5.data[2].y, -0.075 * 2.0, tol);
 248     MY_BOOST_CHECK_SMALL(h_force_5.data[2].z, tol_small);
 249     MY_BOOST_CHECK_CLOSE(h_force_5.data[2].w, 9.375e-4 * 2.0, tol);
 250     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_5.data[0*pitch+2]
 251                                        +h_virial_5.data[3*pitch+2]
 252                                        +h_virial_5.data[5*pitch+2]), -0.02, tol);
 253
 254     MY_BOOST_CHECK_SMALL(h_force_5.data[3].x, tol_small);
 255     MY_BOOST_CHECK_CLOSE(h_force_5.data[3].y, 0.075 * 2.0, tol);
 256     MY_BOOST_CHECK_SMALL(h_force_5.data[3].z, tol_small);
 257     MY_BOOST_CHECK_CLOSE(h_force_5.data[3].w, 9.375e-4 * 2.0, tol);
 258     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_5.data[0*pitch+3]
 259                                        +h_virial_5.data[3*pitch+3]
 260                                        +h_virial_5.data[5*pitch+3]), -0.02, tol);
 261
 262     MY_BOOST_CHECK_SMALL(h_force_5.data[4].x, tol_small);
 263     MY_BOOST_CHECK_SMALL(h_force_5.data[4].y, tol_small);
 264     MY_BOOST_CHECK_CLOSE(h_force_5.data[4].z, -0.45, tol);
 265     MY_BOOST_CHECK_CLOSE(h_force_5.data[4].w, 0.03375, tol);
 266     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_5.data[0*pitch+4]
 267                                        +h_virial_5.data[3*pitch+4]
 268                                        +h_virial_5.data[5*pitch+4]), -0.06, tol);
 269
 270     MY_BOOST_CHECK_SMALL(h_force_5.data[5].x, tol_small);
 271     MY_BOOST_CHECK_SMALL(h_force_5.data[5].y, tol_small);
 272     MY_BOOST_CHECK_CLOSE(h_force_5.data[5].z, 0.45, tol);
 273     MY_BOOST_CHECK_CLOSE(h_force_5.data[5].w, 0.03375, tol);
 274     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_5.data[0*pitch+5]
 275                                        +h_virial_5.data[3*pitch+5]
 276                                        +h_virial_5.data[5*pitch+5]), -0.06, tol);
 277     }
 278
 279     // one more test: this one will test two things:
 280     // 1) That the forces are computed correctly even if the particles are rearranged in memory
 281     // and 2) That two forces can add to the same particle
 282     boost::shared_ptr<SystemDefinition> sysdef_4(new SystemDefinition(4, BoxDim(100.0, 100.0, 100.0), 1, 1, 0, 0, 0, exec_conf));
 283     boost::shared_ptr<ParticleData> pdata_4 = sysdef_4->getParticleData();
 284     pdata_4->setFlags(~PDataFlags(0));
 285
 286     {
 287     ArrayHandle<Scalar4> h_pos(pdata_4->getPositions(), access_location::host, access_mode::readwrite);
 288     ArrayHandle<unsigned int> h_tag(pdata_4->getTags(), access_location::host, access_mode::readwrite);
 289     ArrayHandle<unsigned int> h_rtag(pdata_4->getRTags(), access_location::host, access_mode::readwrite);
 290
 291     // make a square of particles
 292     h_pos.data[0].x = 0.0; h_pos.data[0].y = 0.0; h_pos.data[0].z = 0.0;
 293     h_pos.data[1].x = 1.0; h_pos.data[1].y = 0; h_pos.data[1].z = 0.0;
 294     h_pos.data[2].x = 0; h_pos.data[2].y = 1.0; h_pos.data[2].z = 0.0;
 295     h_pos.data[3].x = 1.0; h_pos.data[3].y = 1.0; h_pos.data[3].z = 0.0;
 296
 297     h_tag.data[0] = 2;
 298     h_tag.data[1] = 3;
 299     h_tag.data[2] = 0;
 300     h_tag.data[3] = 1;
 301     h_rtag.data[h_tag.data[0]] = 0;
 302     h_rtag.data[h_tag.data[1]] = 1;
 303     h_rtag.data[h_tag.data[2]] = 2;
 304     h_rtag.data[h_tag.data[3]] = 3;
 305     }
 306
 307     // build the bond force compute and try it out
 308     boost::shared_ptr<PotentialBondHarmonic> fc_4 = bf_creator(sysdef_4);
 309     fc_4->setParams(0, make_scalar2(1.5, 1.75));
 310     // only add bonds on the left, top, and bottom of the square
 311     sysdef_4->getBondData()->addBondedGroup(Bond(0, 2,3));
 312     sysdef_4->getBondData()->addBondedGroup(Bond(0, 2,0));
 313     sysdef_4->getBondData()->addBondedGroup(Bond(0, 0,1));
 314
 315     fc_4->compute(0);
 316
 317     {
 318     GPUArray<Scalar4>& force_array_6 =  fc_4->getForceArray();
 319     GPUArray<Scalar>& virial_array_6 =  fc_4->getVirialArray();
 320     unsigned int pitch = virial_array_6.getPitch();
 321     ArrayHandle<Scalar4> h_force_6(force_array_6,access_location::host,access_mode::read);
 322     ArrayHandle<Scalar> h_virial_6(virial_array_6,access_location::host,access_mode::read);
 323     // the right two particles shoul only have a force pulling them right
 324     MY_BOOST_CHECK_CLOSE(h_force_6.data[1].x, 1.125, tol);
 325     MY_BOOST_CHECK_SMALL(h_force_6.data[1].y, tol_small);
 326     MY_BOOST_CHECK_SMALL(h_force_6.data[1].z, tol_small);
 327     MY_BOOST_CHECK_CLOSE(h_force_6.data[1].w, 0.2109375, tol);
 328     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_6.data[0*pitch+1]
 329                                        +h_virial_6.data[3*pitch+1]
 330                                        +h_virial_6.data[5*pitch+1]), 0.1875, tol);
 331
 332     MY_BOOST_CHECK_CLOSE(h_force_6.data[3].x, 1.125, tol);
 333     MY_BOOST_CHECK_SMALL(h_force_6.data[3].y, tol_small);
 334     MY_BOOST_CHECK_SMALL(h_force_6.data[3].z, tol_small);
 335     MY_BOOST_CHECK_CLOSE(h_force_6.data[3].w, 0.2109375, tol);
 336     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_6.data[0*pitch+3]
 337                                        +h_virial_6.data[3*pitch+3]
 338                                        +h_virial_6.data[5*pitch+3]), 0.1875, tol);
 339
 340     // the bottom left particle should have a force pulling down and to the left
 341     MY_BOOST_CHECK_CLOSE(h_force_6.data[0].x, -1.125, tol);
 342     MY_BOOST_CHECK_CLOSE(h_force_6.data[0].y, -1.125, tol);
 343     MY_BOOST_CHECK_SMALL(h_force_6.data[0].z, tol_small);
 344     MY_BOOST_CHECK_CLOSE(h_force_6.data[0].w, 0.421875, tol);
 345     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_6.data[0*pitch+0]
 346                                        +h_virial_6.data[3*pitch+0]
 347                                        +h_virial_6.data[5*pitch+0]), 0.375, tol);
 348
 349     // and the top left particle should have a force pulling up and to the left
 350     MY_BOOST_CHECK_CLOSE(h_force_6.data[2].x, -1.125, tol);
 351     MY_BOOST_CHECK_CLOSE(h_force_6.data[2].y, 1.125, tol);
 352     MY_BOOST_CHECK_SMALL(h_force_6.data[2].z, tol_small);
 353     MY_BOOST_CHECK_CLOSE(h_force_6.data[2].w, 0.421875, tol);
 354     MY_BOOST_CHECK_CLOSE(Scalar(1./3.)*(h_virial_6.data[0*pitch+2]
 355                                        +h_virial_6.data[3*pitch+2]
 356                                        +h_virial_6.data[5*pitch+2]), 0.375, tol);
 357     }
 358     }
 359
 360 //! Compares the output of two PotentialBondHarmonics
 361 void bond_force_comparison_tests(bondforce_creator bf_creator1, bondforce_creator bf_creator2, boost::shared_ptr<ExecutionConfiguration> exec_conf)
 362     {
 363     const unsigned int N = 1000;
 364
 365     // create a particle system to sum forces on
 366     // just randomly place particles. We don't really care how huge the bond forces get: this is just a unit test
 367     RandomInitializer rand_init(N, Scalar(0.2), Scalar(0.9), "A");
 368     boost::shared_ptr<SnapshotSystemData> snap = rand_init.getSnapshot();
 369     snap->bond_data.type_mapping.push_back("A");
 370     boost::shared_ptr<SystemDefinition> sysdef(new SystemDefinition(snap, exec_conf));
 371     boost::shared_ptr<ParticleData> pdata = sysdef->getParticleData();
 372     pdata->setFlags(~PDataFlags(0));
 373
 374     boost::shared_ptr<PotentialBondHarmonic> fc1 = bf_creator1(sysdef);
 375     boost::shared_ptr<PotentialBondHarmonic> fc2 = bf_creator2(sysdef);
 376     fc1->setParams(0, make_scalar2(Scalar(300.0), Scalar(1.6)));
 377     fc2->setParams(0, make_scalar2(Scalar(300.0), Scalar(1.6)));
 378
 379     // add bonds
 380     for (unsigned int i = 0; i < N-1; i++)
 381         {
 382         sysdef->getBondData()->addBondedGroup(Bond(0, i, i+1));
 383         }
 384
 385     // compute the forces
 386     fc1->compute(0);
 387     fc2->compute(0);
 388
 389     // verify that the forces are identical (within roundoff errors)
 390     {
 391     GPUArray<Scalar4>& force_array_7 =  fc1->getForceArray();
 392     GPUArray<Scalar>& virial_array_7 =  fc1->getVirialArray();
 393     unsigned int pitch = virial_array_7.getPitch();
 394     ArrayHandle<Scalar4> h_force_7(force_array_7,access_location::host,access_mode::read);
 395     ArrayHandle<Scalar> h_virial_7(virial_array_7,access_location::host,access_mode::read);
 396     GPUArray<Scalar4>& force_array_8 =  fc2->getForceArray();
 397     GPUArray<Scalar>& virial_array_8 =  fc2->getVirialArray();
 398     ArrayHandle<Scalar4> h_force_8(force_array_8,access_location::host,access_mode::read);
 399     ArrayHandle<Scalar> h_virial_8(virial_array_8,access_location::host,access_mode::read);
 400
 401     // compare average deviation between the two computes
 402     double deltaf2 = 0.0;
 403     double deltape2 = 0.0;
 404     double deltav2[6];
 405     for (unsigned int i = 0; i < 6;  i++)
 406         deltav2[i] = 0.0;
 407
 408     for (unsigned int i = 0; i < N; i++)
 409         {
 410         deltaf2 += double(h_force_8.data[i].x - h_force_7.data[i].x) * double(h_force_8.data[i].x - h_force_7.data[i].x);
 411         deltaf2 += double(h_force_8.data[i].y - h_force_7.data[i].y) * double(h_force_8.data[i].y - h_force_7.data[i].y);
 412         deltaf2 += double(h_force_8.data[i].z - h_force_7.data[i].z) * double(h_force_8.data[i].z - h_force_7.data[i].z);
 413         deltape2 += double(h_force_8.data[i].w - h_force_7.data[i].w) * double(h_force_8.data[i].w - h_force_7.data[i].w);
 414         for (unsigned int j = 0; j < 6; j++)
 415             deltav2[j] += double(h_virial_8.data[j*pitch+i] - h_virial_7.data[j*pitch+i]) * double(h_virial_8.data[j*pitch+i] - h_virial_7.data[j*pitch+i]);
 416
 417         // also check that each individual calculation is somewhat close
 418         }
 419     deltaf2 /= double(pdata->getN());
 420     deltape2 /= double(pdata->getN());
 421     for (unsigned int j = 0; j < 6; j++)
 422         deltav2[j] /= double(pdata->getN());
 423     BOOST_CHECK_SMALL(deltaf2, double(tol_small));
 424     BOOST_CHECK_SMALL(deltape2, double(tol_small));
 425     BOOST_CHECK_SMALL(deltav2[0], double(tol_small));
 426     BOOST_CHECK_SMALL(deltav2[1], double(tol_small));
 427     BOOST_CHECK_SMALL(deltav2[2], double(tol_small));
 428     BOOST_CHECK_SMALL(deltav2[3], double(tol_small));
 429     BOOST_CHECK_SMALL(deltav2[4], double(tol_small));
 430     BOOST_CHECK_SMALL(deltav2[5], double(tol_small));
 431     }
 432     }
 433
 434 //! Check ConstForceCompute to see that it operates properly
 435 void const_force_test(boost::shared_ptr<ExecutionConfiguration> exec_conf)
 436     {
 437     // Generate a simple test particle data
 438     boost::shared_ptr<SystemDefinition> sysdef_2(new SystemDefinition(2, BoxDim(1000.0), 1, 0, 0, 0, 0, exec_conf));
 439     boost::shared_ptr<ParticleData> pdata_2 = sysdef_2->getParticleData();
 440     pdata_2->setFlags(~PDataFlags(0));
 441
 442     pdata_2->setPosition(0,make_scalar3(0.0,0.0,0.0));
 443     pdata_2->setPosition(1,make_scalar3(0.9,0.0,0.0));
 444
 445     // Create the ConstForceCompute and check that it works properly
 446     ConstForceCompute fc(sysdef_2, Scalar(-1.3), Scalar(2.5), Scalar(45.67));
 447     {
 448     GPUArray<Scalar4>& force_array_9 =  fc.getForceArray();
 449     GPUArray<Scalar>& virial_array_9 =  fc.getVirialArray();
 450     unsigned int pitch = virial_array_9.getPitch();
 451     ArrayHandle<Scalar4> h_force_9(force_array_9,access_location::host,access_mode::read);
 452     ArrayHandle<Scalar> h_virial_9(virial_array_9,access_location::host,access_mode::read);
 453     MY_BOOST_CHECK_CLOSE(h_force_9.data[0].x, -1.3, tol);
 454     MY_BOOST_CHECK_CLOSE(h_force_9.data[0].y, 2.5, tol);
 455     MY_BOOST_CHECK_CLOSE(h_force_9.data[0].z, 45.67, tol);
 456     MY_BOOST_CHECK_SMALL(h_force_9.data[0].w, tol_small);
 457     MY_BOOST_CHECK_SMALL(h_virial_9.data[0*pitch+0]
 458                         +h_virial_9.data[3*pitch+0]
 459                         +h_virial_9.data[5*pitch+0], tol_small);
 460
 461     MY_BOOST_CHECK_CLOSE(h_force_9.data[1].x, -1.3, tol);
 462     MY_BOOST_CHECK_CLOSE(h_force_9.data[1].y, 2.5, tol);
 463     MY_BOOST_CHECK_CLOSE(h_force_9.data[1].z, 45.67, tol);
 464     MY_BOOST_CHECK_SMALL(h_force_9.data[1].w, tol_small);
 465     MY_BOOST_CHECK_SMALL(h_virial_9.data[0*pitch+1]
 466                         +h_virial_9.data[3*pitch+1]
 467                         +h_virial_9.data[5*pitch+1], tol_small);
 468     }
 469
 470     // check the setforce method
 471     fc.setForce(Scalar(67.54), Scalar(22.1), Scalar(-1.4));
 472     {
 473     GPUArray<Scalar4>& force_array_10 =  fc.getForceArray();
 474     GPUArray<Scalar>& virial_array_10 =  fc.getVirialArray();
 475     unsigned int pitch = virial_array_10.getPitch();
 476     ArrayHandle<Scalar4> h_force_10(force_array_10,access_location::host,access_mode::read);
 477     ArrayHandle<Scalar> h_virial_10(virial_array_10,access_location::host,access_mode::read);
 478     MY_BOOST_CHECK_CLOSE(h_force_10.data[0].x, 67.54, tol);
 479     MY_BOOST_CHECK_CLOSE(h_force_10.data[0].y, 22.1, tol);
 480     MY_BOOST_CHECK_CLOSE(h_force_10.data[0].z, -1.4, tol);
 481     MY_BOOST_CHECK_SMALL(h_force_10.data[0].w, tol_small);
 482     MY_BOOST_CHECK_SMALL(h_virial_10.data[0*pitch+0]
 483                         +h_virial_10.data[3*pitch+0]
 484                         +h_virial_10.data[5*pitch+0], tol_small);
 485
 486     MY_BOOST_CHECK_CLOSE(h_force_10.data[1].x, 67.54, tol);
 487     MY_BOOST_CHECK_CLOSE(h_force_10.data[1].y, 22.1, tol);
 488     MY_BOOST_CHECK_CLOSE(h_force_10.data[1].z, -1.4, tol);
 489     MY_BOOST_CHECK_SMALL(h_force_10.data[1].w, tol_small);
 490     MY_BOOST_CHECK_SMALL(h_virial_10.data[0*pitch+1]
 491                         +h_virial_10.data[3*pitch+1]
 492                         +h_virial_10.data[5*pitch+1], tol_small);
 493     }
 494     }
 495
 496 //! PotentialBondHarmonic creator for bond_force_basic_tests()
 497 boost::shared_ptr<PotentialBondHarmonic> base_class_bf_creator(boost::shared_ptr<SystemDefinition> sysdef)
 498     {
 499     return boost::shared_ptr<PotentialBondHarmonic>(new PotentialBondHarmonic(sysdef));
 500     }
 501
 502 #ifdef ENABLE_CUDA
 503 //! PotentialBondHarmonic creator for bond_force_basic_tests()
 504 boost::shared_ptr<PotentialBondHarmonic> gpu_bf_creator(boost::shared_ptr<SystemDefinition> sysdef)
 505     {
 506     return boost::shared_ptr<PotentialBondHarmonic>(new PotentialBondHarmonicGPU(sysdef));
 507     }
 508 #endif
 509
 510 //! boost test case for bond forces on the CPU
 511 BOOST_AUTO_TEST_CASE( PotentialBondHarmonic_basic )
 512     {
 513     bondforce_creator bf_creator = bind(base_class_bf_creator, _1);
 514     bond_force_basic_tests(bf_creator, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::CPU)));
 515     }
 516
 517 #ifdef ENABLE_CUDA
 518 //! boost test case for bond forces on the GPU
 519 BOOST_AUTO_TEST_CASE( PotentialBondHarmonicGPU_basic )
 520     {
 521     bondforce_creator bf_creator = bind(gpu_bf_creator, _1);
 522     bond_force_basic_tests(bf_creator, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::GPU)));
 523     }
 524
 525 //! boost test case for comparing bond GPU and CPU BondForceComputes
 526 BOOST_AUTO_TEST_CASE( PotentialBondHarmonicGPU_compare )
 527     {
 528     bondforce_creator bf_creator_gpu = bind(gpu_bf_creator, _1);
 529     bondforce_creator bf_creator = bind(base_class_bf_creator, _1);
 530     bond_force_comparison_tests(bf_creator, bf_creator_gpu, boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::GPU)));
 531     }
 532
 533 #endif
 534
 535 //! boost test case for constant forces
 536 BOOST_AUTO_TEST_CASE( ConstForceCompute_basic )
 537     {
 538     const_force_test(boost::shared_ptr<ExecutionConfiguration>(new ExecutionConfiguration(ExecutionConfiguration::CPU)));
 539     }
 540
 541 #ifdef WIN32
 542 #pragma warning( pop )
 543 #endif