         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,

        part1, part2, part3, bothBits);
}

__global__ void densmatr_twoQubitDepolariseKernel(
    REAL depolLevel, REAL* vecReal, REAL *vecImag, long long int numAmpsToVisit,
    long long int part1, long long int part2, long long int part3, 

    vecImag[stateInd] *= fac;
}

// @TODO is separating these 12 amplitudes really faster than letting every 16th base modify 12 elems?
void densmatr_twoQubitDephase(QubitRegister qureg, int qubit1, int qubit2, REAL dephase) {
    if (dephase == 0)
        return;
