         part1, part2, part3, colBit, rowBit);
 }
 
-
-
+/** Called 12 times for every 16 amplitudes in density matrix 
+ * Each sums from the |..0..0..><..0..0..| index to visit either
+ * |..0..0..><..0..1..|,  |..0..0..><..1..0..|,  |..0..0..><..1..1..|,  |..0..1..><..0..0..|
+ * etc and so on to |..1..1..><..1..0|. Labels |part1 0 part2 0 par><t3 0 part4 0 part5|.
+ * From the brain of Simon Benjamin
+ */
 __global__ void densmatr_twoQubitDephaseKernel(
     REAL fac, REAL* vecReal, REAL *vecImag, long long int numBackgroundStates, long long int numAmpsToVisit,
     long long int part1, long long int part2, long long int part3, long long int part4, long long int part5,

    CUDABlocks = ceil((REAL)(qureg.numAmpsPerChunk>>1)/threadsPerCUDABlock);
    statevec_sigmaXKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, targetQubit);
}


__global__ void statevec_sigmaYKernel(QubitRegister qureg, const int targetQubit, const int conjFac){


    // initialise the state to all zeros
    long long int index = blockIdx.x*blockDim.x + threadIdx.x;
    if (index >= densityNumElems) return;
    
    densityReal[index] = 0.0;
    densityImag[index] = 0.0;
    
