     cudaMemcpyToSymbol(GPU_rho_coeff, &(CPU_rho_coeff[0]), order * (2*order+1) * sizeof(float));
 
     // setup the grid to run the kernel with one thread per particle in the group
+    dim3 grid( (int)ceil((double)group_size / (double)block_size), 1, 1);
+    dim3 threads(block_size, 1, 1);
+
+    // setup the grid to run the kernel with one thread per particle in the group
     dim3 P_grid( (int)ceil((double)group_size / (double)block_size), 1, 1);
     dim3 P_threads(block_size, 1, 1);
     

        local_field.z = E_z[tid].x;
      
        Electric_field[tid] = local_field;
        }
    }


    if (error != cudaSuccess)
        return error;
        
    // set the grid charge to zero
    cudaMemset(GPU_rho_real_space, 0.0f, sizeof(cufftComplex)*Nx*Ny*Nz);

    // zero the force arrays for all particles
//    zero_forces <<< grid, threads >>> (force_data, pdata);
    cudaMemset(force_data.force, 0.0f, sizeof(float4)*pdata.N);
    cudaMemset(force_data.virial, 0.0f, sizeof(float)*pdata.N);


    // run the kernels
    // assign charges to the grid points, one thread per particles
