     int numThreadsPerBlock = 128;
     int numBlocks = ceil(qureg.numAmpsPerChunk / (qreal) numThreadsPerBlock);
     
-    
     // create global GPU array for outcomeProbs
     qreal* d_outcomeProbs;
     long long int numOutcomes = (1LL << numQubits);

    cudaMemcpy(d_termCoeffs, hamil.termCoeffs, mem_termCoeffs, cudaMemcpyHostToDevice);
    
    int numThreadsPerBlock = 128;
    int numBlocks = ceil(op.numElemsPerChunk)/numThreadsPerBlock);
    agnostic_initDiagonalOpFromPauliHamilKernel<<numBlocks, numThreadsPerBlock>>(
        op, d_pauliCodes, d_termCoeffs, hamil.numSumTerms);
    
    // copy populated operator into to RAM
    cudaDeviceSynchronize();
    cudaMemcpy(op.real, op.deviceOperator.real, op.numElemsPerChunk, cudaMemcpyDeviceToHost);
    cudaMemcpy(op.imag, op.deviceOperator.imag, op.numElemsPerChunk, cudaMemcpyDeviceToHost);

    cudaFree(d_pauliCodes);
    cudaFree(d_termCoeffs);


void agnostic_syncDiagonalOp(DiagonalOp op) {
    cudaDeviceSynchronize();
    size_t mem_elems = op.numElemsPerChunk * sizeof *op.real;
    cudaMemcpy(op.deviceOperator.real, op.real, mem_elems, cudaMemcpyHostToDevice);
    cudaMemcpy(op.deviceOperator.imag, op.imag, mem_elems, cudaMemcpyHostToDevice);
}

__global__ void agnostic_initDiagonalOpFromPauliHamilKernel(
