 #include "hoomd/TextureTools.h"
 #include "hoomd/GPUPartition.cuh"
 
-#ifdef WIN32
-#include <cassert>
-#else
-#include <assert.h>
-#endif
+#ifdef NVCC
+#include "hoomd/WarpTools.cuh"
+#endif // NVCC
 
 /*! \file AnisoPotentialPairGPU.cuh
     \brief Defines templated GPU kernel code for calculating the anisotropic ptl pair forces and torques

                {
                case 0:
                    {
                    static unsigned int max_block_size = UINT_MAX;
                    if (max_block_size == UINT_MAX)
                        max_block_size = aniso_get_max_block_size(gpu_compute_pair_aniso_forces_kernel<evaluator, 0, 0>);


                    block_size = block_size < max_block_size ? block_size : max_block_size;
                    dim3 grid(nwork / (block_size/tpp) + 1, 1, 1);

                    shared_bytes += sizeof(Scalar)*block_size;

                    gpu_compute_pair_aniso_forces_kernel<evaluator, 0, 0>
                        <<<grid, block_size, shared_bytes>>>(pair_args.d_force,
                                                      pair_args.d_torque,
                                                      pair_args.d_virial,
                                                      pair_args.virial_pitch,
                                                      nwork,
                                                      pair_args.d_pos,
                                                      pair_args.d_diameter,
                                                      pair_args.d_charge,
                                                      pair_args.d_orientation,
                                                      pair_args.box,
                                                      pair_args.d_n_neigh,
                                                      pair_args.d_nlist,
                                                      pair_args.d_head_list,
                                                      d_params,
                                                      pair_args.d_rcutsq,
                                                      pair_args.ntypes,
                                                      tpp,
                                                      range.first);
                    break;
                    }
                case 1:
                    {
                    static unsigned int max_block_size = UINT_MAX;
                    if (max_block_size == UINT_MAX)
                        max_block_size = aniso_get_max_block_size(gpu_compute_pair_aniso_forces_kernel<evaluator, 1, 0>);

                    block_size = block_size < max_block_size ? block_size : max_block_size;
                    dim3 grid(nwork / (block_size/tpp) + 1, 1, 1);

                    shared_bytes += sizeof(Scalar)*block_size;

                    gpu_compute_pair_aniso_forces_kernel<evaluator, 1, 0>
                        <<<grid, block_size, shared_bytes>>>(pair_args.d_force,
                                                      pair_args.d_torque,
                                                      pair_args.d_virial,
                                                      pair_args.virial_pitch,
                                                      nwork,
                                                      pair_args.d_pos,
                                                      pair_args.d_diameter,
                                                      pair_args.d_charge,
                                                      pair_args.d_orientation,
                                                      pair_args.box,
                                                      pair_args.d_n_neigh,
                                                      pair_args.d_nlist,
                                                      pair_args.d_head_list,
                                                      d_params,
                                                      pair_args.d_rcutsq,
                                                      pair_args.ntypes,
                                                      tpp,
                                                      range.first);
                    break;
                    }
                default:

    Scalar virialyz = 0.0f;
    Scalar virialzz = 0.0f;

    if (active)
        {
        // load in the length of the neighbor list
        unsigned int n_neigh = d_n_neigh[idx];

        // read in the position of our particle
        Scalar4 postypei = __ldg(d_pos + idx);
        Scalar3 posi = make_scalar3(postypei.x, postypei.y, postypei.z);
        Scalar4 quati = __ldg(d_orientation + idx);

        Scalar di = 0.0f;
        if (evaluator::needsDiameter())
            di = __ldg(d_diameter + idx);

        Scalar qi = 0.0f;
        if (evaluator::needsCharge())
            qi = __ldg(d_charge + idx);

        unsigned int my_head = d_head_list[idx];
        unsigned int cur_j = 0;

        unsigned int next_j(0);
        next_j = threadIdx.x%tpp < n_neigh ? __ldg(d_nlist + my_head + threadIdx.x%tpp) : 0;

        // loop over neighbors
        for (int neigh_idx = threadIdx.x%tpp; neigh_idx < n_neigh; neigh_idx+=tpp)
            {
                {
                // read the current neighbor index
                // prefetch the next value and set the current one
                cur_j = next_j;
                if (neigh_idx+tpp < n_neigh)
                    {
                    next_j = __ldg(d_nlist + my_head + neigh_idx+tpp);
                    }

                // get the neighbor's position
                Scalar4 postypej = __ldg(d_pos + cur_j);
                Scalar3 posj = make_scalar3(postypej.x, postypej.y, postypej.z);
                Scalar4 quatj = __ldg(d_orientation + cur_j);

                Scalar dj = 0.0f;
                if (evaluator::needsDiameter())
                    dj = __ldg(d_diameter + cur_j);

                Scalar qj = 0.0f;
                if (evaluator::needsCharge())
                    qj = __ldg(d_charge + cur_j);

                // calculate dr (with periodic boundary conditions)
                Scalar3 dx = posi - posj;

                // apply periodic boundary conditions
                dx = box.minImage(dx);

                // calculate r squared
                Scalar rsq = dot(dx, dx);

                // access the per type pair parameters
                unsigned int typpair = typpair_idx(__scalar_as_int(postypei.w), __scalar_as_int(postypej.w));
                Scalar rcutsq = s_rcutsq[typpair];
                typename evaluator::param_type param = s_params[typpair];

                // design specifies that energies are shifted if
                // 1) shift mode is set to shift
                bool energy_shift = false;
                if (shift_mode == 1)
                    energy_shift = true;

                // evaluate the potential
                Scalar3 jforce = { 0.0f, 0.0f, 0.0f };
                Scalar3 torquei = { 0.0f, 0.0f, 0.0f };
                Scalar3 torquej = { 0.0f, 0.0f, 0.0f };
                Scalar pair_eng = 0.0f;

                // constructor call
                evaluator eval(dx, quati, quatj, rcutsq, param);
                if (evaluator::needsDiameter())
                    eval.setDiameter(di, dj);
                if (evaluator::needsCharge())
                    eval.setCharge(qi, qj);

                // call evaluator
                eval.evaluate(jforce, pair_eng, energy_shift, torquei, torquej);

                // calculate the virial
                if (compute_virial)
                    {
                    Scalar3 jforce2 = Scalar(0.5)*jforce;
                    virialxx +=  dx.x * jforce2.x;
                    virialxy +=  dx.y * jforce2.x;
                    virialxz +=  dx.z * jforce2.x;
                    virialyy +=  dx.y * jforce2.y;
                    virialyz +=  dx.z * jforce2.y;
                    virialzz +=  dx.z * jforce2.z;
                    }

                // add up the force vector components
                force.x += jforce.x;
                force.y += jforce.y;
                force.z += jforce.z;
                torque.x += torquei.x;
                torque.y += torquei.y;
                torque.z += torquei.z;

                force.w += pair_eng;
                }
            }

        // potential energy per particle must be halved
        force.w *= Scalar(0.5);
        }

    // reduce force over threads in cta
    hoomd::detail::WarpReduce<Scalar, tpp> reducer;
    force.x = reducer.Sum(force.x);
    force.y = reducer.Sum(force.y);
    force.z = reducer.Sum(force.z);
    force.w = reducer.Sum(force.w);

    torque.x = reducer.Sum(torque.x);
    torque.y = reducer.Sum(torque.y);
    torque.z = reducer.Sum(torque.z);

    // now that the force calculation is complete, write out the result
    if (active && threadIdx.x % tpp == 0)
        {
        d_force[idx] = force;
        d_torque[idx] = torque;
