         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
     );
 }
 

//==                            (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
(

}


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