 //==                            (cuDoubleComplex*)beta_phi);
 //== }
 
-//== __global__ void restore_valence_density_gpu_kernel(int num_gvec_loc,
-//==                                                    int* atom_type,
-//==                                                    int* num_beta, 
-//==                                                    double* atom_pos,
-//==                                                    int* gvec,
-//==                                                    cuDoubleComplex* pp_complex_density_matrix,
-//==                                                    int ldm,
-//==                                                    cuDoubleComplex** q_pw,
-//==                                                    cuDoubleComplex* f_pw)
-//== {
-//==     extern __shared__ char sdata_ptr[];
-//==     cuDoubleComplex* sdata = (cuDoubleComplex*)&sdata_ptr[0];
-//== 
-//==     int ia = blockIdx.x;
-//== 
-//==     int iat = atom_type[ia];
-//== 
-//==     int nbf = num_beta[iat];
-//== 
-//==     cuDoubleComplex* q_pw_t = q_pw[iat];
-//==     //printf("ia : %i, type : %i, nbf : %i, q_pw : %p", ia, iat, nbf, q_pw_t);
-//== 
-//==     double ax = atom_pos[array2D_offset(0, ia, 3)];
-//==     double ay = atom_pos[array2D_offset(1, ia, 3)];
-//==     double az = atom_pos[array2D_offset(2, ia, 3)];
-//== 
-//==     if (threadIdx.x == 0)
-//==     {
-//==         for (int xi2 = 0; xi2 < nbf; xi2++)
-//==         {
-//==             for (int xi1 = 0; xi1 <= xi2; xi1++)
-//==             {
-//==                 int idx12 = xi2 * (xi2 + 1) / 2 + xi1;
-//==                 sdata[idx12] = pp_complex_density_matrix[array4D_offset(xi2, xi1, 0, ia, ldm, ldm, 1)];
-//==             }
-//==         }
-//==     }
-//==     __syncthreads();
-//== 
-//==     cuDoubleComplex* f_pw_a = &f_pw[array2D_offset(0, ia, num_gvec_loc)];
-//==     
-//==     int N = num_blocks(num_gvec_loc, blockDim.x);
-//== 
-//==     for (int n = 0; n < N; n++)
-//==     {
-//==         int igloc = n * blockDim.x + threadIdx.x;
-//==         if (igloc < num_gvec_loc)
-//==         {
-//==             int gvx = gvec[array2D_offset(0, igloc, 3)];
-//==             int gvy = gvec[array2D_offset(1, igloc, 3)];
-//==             int gvz = gvec[array2D_offset(2, igloc, 3)];
-//== 
-//==             double p = twopi * (ax * gvx + ay * gvy + az * gvz);
-//==             
-//==             double sinp = sin(p);
-//==             double cosp = cos(p);
-//== 
-//==             cuDoubleComplex zval = make_cuDoubleComplex(0.0, 0.0);
-//== 
-//==             // \sum_{xi1, xi2} D_{xi2,xi1} * Q(G)_{xi1, xi2}
-//==             for (int xi2 = 0; xi2 < nbf; xi2++)
-//==             {
-//==                 int idx12 = xi2 * (xi2 + 1) / 2;
-//== 
-//==                 //cuDoubleComplex q = cuCmul(make_cuDoubleComplex(cosp, -sinp), q_pw_t[array2D_offset(igloc, idx12 + xi2, num_gvec_loc)]);
-//== 
-//==                 // add diagonal term
-//==                 //f_pw_a[igloc] = cuCadd(f_pw_a[igloc], cuCmul(sdata[idx12 + xi2], q));
-//==                 zval = cuCadd(zval, cuCmul(sdata[idx12 + xi2], q_pw_t[array2D_offset(igloc, idx12 + xi2, num_gvec_loc)]));
-//== 
-//==                 // add non-diagonal terms
-//==                 for (int xi1 = 0; xi1 < xi2; xi1++, idx12++)
-//==                 {
-//==                     cuDoubleComplex q = q_pw_t[array2D_offset(igloc, idx12, num_gvec_loc)];
-//==                     //q = cuCmul(make_cuDoubleComplex(cosp, -sinp), q_pw_t[array2D_offset(igloc, idx12, num_gvec_loc)]);
-//==                     
-//==                     //double d = 2 * cuCreal(cuCmul(sdata[idx12], q));
-//== 
-//==                     //f_pw_a[igloc] = cuCadd(f_pw_a[igloc], make_cuDoubleComplex(d, 0));
-//==                     //double d = 2 * cuCreal(cuCmul(sdata[idx12], q_pw_t[array2D_offset(igloc, idx12, num_gvec_loc)])
-//==                     zval.x += 2 * (sdata[idx12].x * q.x - sdata[idx12].y * q.y);
-//==                     //zval = cuCadd(zval, make_cuDoubleComplex(2 * cuCreal(cuCmul(sdata[idx12], q_pw_t[array2D_offset(igloc, idx12, num_gvec_loc)])), 0.0));
-//==                 }
-//==             }
-//==             f_pw_a[igloc] = cuCadd(f_pw_a[igloc], cuCmul(zval, make_cuDoubleComplex(cosp, -sinp))); 
-//==         }
-//==     }
-//== }
-//== 
-//== __global__ void reduce_rho_pw_kernel(int num_atoms, int num_gvec_loc, cuDoubleComplex* f_pw, cuDoubleComplex* rho_pw)
-//== {
-//==     int igloc = blockDim.x * blockIdx.x + threadIdx.x;
-//== 
-//==     if (igloc < num_gvec_loc)
-//==     {
-//==         for (int ia = 0; ia < num_atoms; ia++) 
-//==             rho_pw[igloc] = cuCadd(rho_pw[igloc], f_pw[array2D_offset(igloc, ia, num_gvec_loc)]);
-//==     }
-//== }
-//== 
-//== 
-//== extern "C" void restore_valence_density_gpu(int num_atoms, 
-//==                                             int num_gvec_loc,
-//==                                             int* atom_type,
-//==                                             int* num_beta, 
-//==                                             double* atom_pos, 
-//==                                             int* gvec,
-//==                                             void* pp_complex_density_matrix,
-//==                                             int ldm,
-//==                                             void** q_pw,
-//==                                             void* rho_pw)
-//== {
-//==     dim3 threadsPerBlock(1024);
-//==     dim3 numBlocks(num_atoms);
-//== 
-//==     cuDoubleComplex* f_pw;
-//==     f_pw = (cuDoubleComplex*)cuda_malloc(num_gvec_loc * num_atoms * sizeof(cuDoubleComplex));
-//==     cuda_memset(f_pw, 0, num_gvec_loc * num_atoms * sizeof(cuDoubleComplex));
-//== 
-//==     restore_valence_density_gpu_kernel<<<
-//==         numBlocks,
-//==         threadsPerBlock,
-//==         sizeof(cuDoubleComplex) * ldm * (ldm + 1) / 2>>>(num_gvec_loc,
-//==                                                          atom_type,
-//==                                                          num_beta, 
-//==                                                          atom_pos, 
-//==                                                          gvec, 
-//==                                                          (cuDoubleComplex*)pp_complex_density_matrix,
-//==                                                          ldm,
-//==                                                          (cuDoubleComplex**)q_pw,
-//==                                                          f_pw);
-//==     
-//==     cuda_memset(rho_pw, 0, num_gvec_loc * sizeof(cuDoubleComplex));
-//==     
-//==     dim3 grid_t(128);
-//==     dim3 grid_b(num_blocks(num_gvec_loc, grid_t.x));
-//==     reduce_rho_pw_kernel<<<grid_b, grid_t>>>
-//==         (num_atoms, num_gvec_loc, f_pw, (cuDoubleComplex*)rho_pw);
-//==     
-//==     cuda_device_synchronize();
-//==     cuda_free(f_pw);
-//== }
-
-
-
-
-//== __global__ void restore_valence_density_gpu_kernel_v2
-//== (
-//==     int num_gvec_loc,
-//==     int num_beta, 
-//==     double ax,
-//==     double ay,
-//==     double az,
-//==     int* gvec,
-//==     cuDoubleComplex* pp_complex_density_matrix,
-//==     int ldm,
-//==     cuDoubleComplex* q_pw_t,
-//==     cuDoubleComplex* rho_pw
-//== )
-//== {
-//==     extern __shared__ char sdata_ptr[];
-//==     cuDoubleComplex* sdata = (cuDoubleComplex*)&sdata_ptr[0];
-//== 
-//==     if (threadIdx.x == 0)
-//==     {
-//==         for (int xi2 = 0; xi2 < num_beta; xi2++)
-//==         {
-//==             for (int xi1 = 0; xi1 <= xi2; xi1++)
-//==             {
-//==                 int idx12 = xi2 * (xi2 + 1) / 2 + xi1;
-//==                 sdata[idx12] = pp_complex_density_matrix[array3D_offset(xi2, xi1, 0, ldm, ldm)];
-//==             }
-//==         }
-//==     }
-//==     __syncthreads();
-//== 
-//==     int igloc = blockIdx.x * blockDim.x + threadIdx.x;
-//==     if (igloc < num_gvec_loc)
-//==     {
-//==         int gvx = gvec[array2D_offset(0, igloc, 3)];
-//==         int gvy = gvec[array2D_offset(1, igloc, 3)];
-//==         int gvz = gvec[array2D_offset(2, igloc, 3)];
-//== 
-//==         double p = twopi * (ax * gvx + ay * gvy + az * gvz);
-//==         
-//==         double sinp = sin(p);
-//==         double cosp = cos(p);
-//== 
-//==         cuDoubleComplex zval = make_cuDoubleComplex(0.0, 0.0);
-//== 
-//==         // \sum_{xi1, xi2} D_{xi2,xi1} * Q(G)_{xi1, xi2}
-//==         for (int xi2 = 0; xi2 < num_beta; xi2++)
-//==         {
-//==             int idx12 = xi2 * (xi2 + 1) / 2;
-//== 
-//==             // add diagonal term
-//==             zval = cuCadd(zval, cuCmul(sdata[idx12 + xi2], q_pw_t[array2D_offset(igloc, idx12 + xi2, num_gvec_loc)]));
-//== 
-//==             // add non-diagonal terms
-//==             for (int xi1 = 0; xi1 < xi2; xi1++, idx12++)
-//==             {
-//==                 cuDoubleComplex q = q_pw_t[array2D_offset(igloc, idx12, num_gvec_loc)];
-//==                 zval.x += 2 * (sdata[idx12].x * q.x - sdata[idx12].y * q.y);
-//==             }
-//==         }
-//==         rho_pw[igloc] = cuCadd(rho_pw[igloc], cuCmul(zval, make_cuDoubleComplex(cosp, -sinp))); 
-//==     }
-//== }
-//== 
-//== extern "C" void restore_valence_density_gpu_v2(int num_gvec_loc,
-//==                                                int num_beta,
-//==                                                double ax,
-//==                                                double ay,
-//==                                                double az,
-//==                                                int* gvec,
-//==                                                void* pp_complex_density_matrix,
-//==                                                int ldm,
-//==                                                void* q_pw_t,
-//==                                                void* rho_pw,
-//==                                                int stream_id)
-//== {
-//==     cudaStream_t stream = (stream_id == -1) ? NULL : streams[stream_id];
-//== 
-//==     dim3 grid_t(64);
-//==     dim3 grid_b(num_blocks(num_gvec_loc, grid_t.x));
-//== 
-//==     restore_valence_density_gpu_kernel_v2<<<grid_b, grid_t, sizeof(cuDoubleComplex) * ldm * (ldm + 1) / 2, stream>>>
-//==         (num_gvec_loc, num_beta, ax, ay, az, gvec, (cuDoubleComplex*)pp_complex_density_matrix, ldm,
-//==          (cuDoubleComplex*)q_pw_t, (cuDoubleComplex*)rho_pw);
-//== }
-
-//==__global__ void compute_d_mtrx_gpu_kernel
-//==(
-//==    int num_gvec_loc, 
-//==    cuDoubleComplex const* vtmp, 
-//==    cuDoubleComplex const* q_pw_t, 
-//==    cuDoubleComplex* d_mtrx_gpu
-//==)
-//=={
-//==    int idx = blockIdx.x;
-//==
-//==    //int N = num_blocks(num_gvec_loc, blockDim.x);
-//==
-//==    extern __shared__ char sdata_ptr[];
-//==    cuDoubleComplex* sdata = (cuDoubleComplex*)&sdata_ptr[0];
-//==
-//==    //sdata[threadIdx.x] = make_cuDoubleComplex(0.0, 0.0);
-//==
-//==    cuDoubleComplex z = make_cuDoubleComplex(0, 0);
-//==    int igloc = threadIdx.x;
-//==    while (igloc < num_gvec_loc)
-//==    {
-//==        z = cuCadd(z, cuCmul(vtmp[igloc], cuConj(q_pw_t[array2D_offset(igloc, idx,  num_gvec_loc)])));
-//==        igloc += blockDim.x;
-//==    }
-//==    sdata[threadIdx.x] = z;
-//==    __syncthreads();
-//==
-//==    //== for (int n = 0; n < N; n++)
-//==    //== {
-//==    //==     int igloc = n * blockDim.x + threadIdx.x;
-//==    //==     if (igloc < num_gvec_loc)
-//==    //==     {
-//==    //==         sdata[threadIdx.x] = cuCadd(sdata[threadIdx.x], 
-//==    //==                                     cuCmul(vtmp[igloc], 
-//==    //==                                            cuConj(q_pw_t[array2D_offset(igloc, idx,  num_gvec_loc)])));
-//==    //==     }
-//==    //== }
-//==    //== 
-//==    //== __syncthreads();
-//==
-//==    for (int s = 1; s < blockDim.x; s *= 2) 
-//==    {
-//==        if (threadIdx.x % (2 * s) == 0) sdata[threadIdx.x] = cuCadd(sdata[threadIdx.x], sdata[threadIdx.x + s]);
-//==        __syncthreads();
-//==    }
-//==
-//==    d_mtrx_gpu[idx] = sdata[0];
-//==}
-
-//== extern "C" void compute_d_mtrx_valence_gpu(int num_gvec_loc,
-//==                                            int num_elements,
-//==                                            cuDoubleComplex const* veff, 
-//==                                            int const* gvec, 
-//==                                            double ax,
-//==                                            double ay,
-//==                                            double az,
-//==                                            cuDoubleComplex* vtmp,
-//==                                            cuDoubleComplex const* q_pw_t,
-//==                                            cuDoubleComplex* d_mtrx,
-//==                                            int stream_id)
-//== {
-//==     cudaStream_t stream = (stream_id == -1) ? NULL : streams[stream_id];
-//== 
-//==     dim3 grid_t(128);
-//== 
-//==     dim3 grid_b(num_blocks(num_gvec_loc, grid_t.x));
-//== 
-//==     mul_veff_with_phase_factors_gpu_kernel <<<grid_b, grid_t, 0, stream>>>
-//==     (
-//==         num_gvec_loc,
-//==         veff,
-//==         gvec,
-//==         ax,
-//==         ay,
-//==         az, 
-//==         vtmp
-//==     );
-//== 
-//==     grid_b = dim3(num_elements);
-//== 
-//==     compute_d_mtrx_gpu_kernel <<<grid_b, grid_t, grid_t.x * sizeof(cuDoubleComplex), stream>>>
-//==     (
-//==         num_gvec_loc,
-//==         vtmp,
-//==         q_pw_t,
-//==         d_mtrx
-//==     );
-//== 
-//== }
-
-
-
-// TODO: proper order of indices for atom_pos
-
-__global__ void generate_phase_factors_conj_gpu_kernel
-(
-    int num_gvec_loc, 
-    int num_atoms, 
-    double* atom_pos, 
-    int* gvec, 
-    cuDoubleComplex* phase_factors
-)
-{
-    int ia = blockIdx.y;
-    int igloc = blockIdx.x * blockDim.x + threadIdx.x;
-
-    if (igloc < num_gvec_loc)
-    {
-        int gvx = gvec[array2D_offset(0, igloc, 3)];
-        int gvy = gvec[array2D_offset(1, igloc, 3)];
-        int gvz = gvec[array2D_offset(2, igloc, 3)];
-    
-        double ax = atom_pos[array2D_offset(ia, 0, num_atoms)];
-        double ay = atom_pos[array2D_offset(ia, 1, num_atoms)];
-        double az = atom_pos[array2D_offset(ia, 2, num_atoms)];
-
-        double p = twopi * (ax * gvx + ay * gvy + az * gvz);
-
-        double sinp = sin(p);
-        double cosp = cos(p);
-
-        phase_factors[array2D_offset(igloc, ia, num_gvec_loc)] = make_cuDoubleComplex(cosp, -sinp);
-    }
-}
+
 
 __global__ void generate_phase_factors_gpu_kernel
 (

        atom_pos__, 
        gvec__, 
        phase_factors__
    );
}

__global__ void sum_q_pw_d_mtrx_pw_gpu_kernel
(
    int num_gvec_loc,
    int num_beta,
    cuDoubleComplex* q_pw_t,
    cuDoubleComplex* d_mtrx_pw,
    cuDoubleComplex* rho_pw
)
{
    int igloc = blockIdx.x * blockDim.x + threadIdx.x;
    if (igloc < num_gvec_loc)
    {
        cuDoubleComplex zval = make_cuDoubleComplex(0.0, 0.0);

        // \sum_{xi1, xi2} D_{xi2,xi1} * Q(G)_{xi1, xi2}
        for (int xi2 = 0; xi2 < num_beta; xi2++)
        {
            int idx12 = xi2 * (xi2 + 1) / 2;

            // add diagonal term
            zval = cuCadd(zval, cuCmul(d_mtrx_pw[array2D_offset(igloc, xi2 * num_beta + xi2, num_gvec_loc)], 
                                       q_pw_t[array2D_offset(igloc, idx12 + xi2, num_gvec_loc)]));

            // add non-diagonal terms
            for (int xi1 = 0; xi1 < xi2; xi1++, idx12++)
            {
                cuDoubleComplex q = q_pw_t[array2D_offset(igloc, idx12, num_gvec_loc)];
                cuDoubleComplex d1 = d_mtrx_pw[array2D_offset(igloc, xi2 * num_beta + xi1, num_gvec_loc)];
                cuDoubleComplex d2 = d_mtrx_pw[array2D_offset(igloc, xi1 * num_beta + xi2, num_gvec_loc)];

                zval = cuCadd(zval, cuCmul(q, d1));
                zval = cuCadd(zval, cuCmul(cuConj(q), d2));
            }
        }
        rho_pw[igloc] = cuCadd(rho_pw[igloc], zval);
    }
}

extern "C" void sum_q_pw_d_mtrx_pw_gpu(int num_gvec_loc,
                                       int num_beta,
                                       cuDoubleComplex* q_pw_t,
                                       cuDoubleComplex* d_mtrx_pw,
                                       cuDoubleComplex* rho_pw)
{
    CUDA_timer t("sum_q_pw_d_mtrx_pw_gpu");

    dim3 grid_t(64);
    dim3 grid_b(num_blocks(num_gvec_loc, grid_t.x));
    
    sum_q_pw_d_mtrx_pw_gpu_kernel <<<grid_b, grid_t>>>
    (
        num_gvec_loc, 
        num_beta, 
        q_pw_t, 
        d_mtrx_pw, 
        rho_pw
    );
}


}


extern "C" void generate_phase_factors_gpu(int num_gvec_loc__,
                                           int num_atoms__,
                                           int const* gvec__,
