     int ld1 = max_num_mt_points * 4;
     int ld2 = ld1 * (lmax_pw + 1);
     int ld3 = ld2 * num_atom_types;
-    double* bessel_ptr = &bessel[l * ld1 + iat * ld2 + igk * ld3];
-    double* lo_ptr = &lo[ld1 * ilo];
+    double* jl_ptr = &jl_coefs[l * ld1 + iat * ld2 + igk * ld3];
+    double* lo_ptr = &lo_coefs[ld1 * ilo];
     double* r_dr_ptr = &r_dr[2 * max_num_mt_points * iat];
     
-    result[igk + ilo * num_gkvec] = spline_inner_product_gpu_function(max_num_mt_points, nmtp, r_dr_ptr, bessel_ptr, lo_ptr);
-}
-
-
-void bessel_lo_inner_product_gpu(int num_gkvec, int num_rowcol, int max_num_mt_points, int lmax_pw, int num_atom_types, 
-                                 double* bessel, double* lo, int* l_by_ilo, int* iat_by_ilo, int* nmtp_by_ilo, 
-                                 double* r_dr, double* result)
+    jlo[igk + ilo * num_gkvec] = spline_inner_product_gpu_function(max_num_mt_points, nmtp, r_dr_ptr, jl_ptr, lo_ptr);
+}
+
+
+void bessel_lo_inner_product_gpu(int num_gkvec, int num_lo, int max_num_mt_points, int lmax_pw, int num_atom_types, 
+                                 double* jl_coefs, double* lo_coefs, int* l_by_ilo, int* iat_by_ilo, int* nmtp_by_ilo, 
+                                 double* r_dr, double* jlo)
 {
     dim3 threadsPerBlock(64);
-    dim3 numBlocks(num_gkvec, num_rowcol);
+    dim3 numBlocks(num_gkvec, num_lo);
 
     bessel_lo_inner_product_gpu_kernel<<<numBlocks, threadsPerBlock, 64 * 16>>>
-        (max_num_mt_points, lmax_pw, num_atom_types, num_gkvec, bessel, lo, l_by_ilo, iat_by_ilo, nmtp_by_ilo,
-         r_dr, result);
-
-
+        (max_num_mt_points, lmax_pw, num_atom_types, num_gkvec, jl_coefs, lo_coefs, l_by_ilo, iat_by_ilo, nmtp_by_ilo,
+         r_dr, jlo);
+}
+
+__global__ void bessel_vlo_inner_product_gpu_kernel(int max_num_mt_points, int lmax_pw, int lmmax_pw, int num_atom_types, 
+                                                    int num_gkvec, double* jl_coefs, cuDoubleComplex* vlo_coefs,
+                                                    int* l_by_lm, int* iat_by_ilo, int* nmtp_by_ilo, double* r_dr, 
+                                                    cuDoubleComplex* jvlo)
+
+{
+    int igk = blockIdx.x;
+    int ilo = blockIdx.y;
+    int lm = blockIdx.z;
+    int l = l_by_lm[lm];
+    int iat = iat_by_ilo[ilo];
+    int nmtp = nmtp_by_ilo[ilo];
+    
+    int ld1 = max_num_mt_points * 4;
+    int ld2 = ld1 * (lmax_pw + 1);
+    int ld3 = ld2 * num_atom_types;
+    double* jl_ptr = &jl_coefs[l * ld1 + iat * ld2 + igk * ld3];
+    cuDoubleComplex* vlo_ptr = &vlo_coefs[ld1 * lm + ld1 * lmmax_pw * ilo];
+    double* r_dr_ptr = &r_dr[2 * max_num_mt_points * iat];
+    
+    jvlo[lm + igk * lmmax_pw + ilo * lmmax_pw * num_gkvec] = 
+        spline_inner_product_gpu_function(max_num_mt_points, nmtp, r_dr_ptr, jl_ptr, vlo_ptr);
+}
+
+
+
+void bessel_vlo_inner_product_gpu(int num_gkvec, int num_lo, int max_num_mt_points, int lmax_pw, int lmmax_pw, int num_atom_types, 
+                                 double* jl_coefs, void* vlo_coefs, int* l_by_lm, int* iat_by_ilo, int* nmtp_by_ilo, 
+                                 double* r_dr, void* jvlo)
+{
+    dim3 threadsPerBlock(128);
+    dim3 numBlocks(num_gkvec, num_lo, lmmax_pw);
+
+    bessel_vlo_inner_product_gpu_kernel<<<numBlocks, threadsPerBlock, 128 * 16>>>
+        (max_num_mt_points, lmax_pw, lmmax_pw, num_atom_types, num_gkvec, jl_coefs, (cuDoubleComplex*)vlo_coefs, 
+         l_by_lm, iat_by_ilo, nmtp_by_ilo, r_dr, (cuDoubleComplex*)jvlo);
 }
 
 __global__ void bessel_vlm_inner_product_gpu_kernel(int max_num_mt_points, int lmax_pot, int lmmax_pot, int* iat_by_ia,

    return sdata[0];
}

template <typename T>
__global__ void spline_inner_product_gpu_kernel(int ld, int size, double* r_dr, T* s1_coefs, T* s2_coefs, T* result)
{
    result[0] = spline_inner_product_gpu_function(ld, size, r_dr, s1_coefs, s2_coefs);
}


    cudaFree(d_result);
    free(h_result);
    
    //cuDoubleComplex* d_zresult;
    //cudaMalloc(&d_zresult, 1 * sizeof(cuDoubleComplex));
    //
    //cuDoubleComplex* zs2;
    //cudaMalloc(&zs2, size * 4 * sizeof(cuDoubleComplex));
    //
    //for (int i = 0; i < size * 4; i++) zs2[i] = make_cuDoubleComplex(s2_coefs[i], s2_coefs[i]);

    //spline_inner_product_gpu_kernel<<<numBlocks, threadsPerBlock, 64 * 16>>>(size, size, r_dr, s1_coefs, zs2, d_zresult);

    //cuDoubleComplex* h_zresult = (cuDoubleComplex*)malloc(1 * sizeof(cuDoubleComplex));
    //cudaMemcpy(h_zresult, d_zresult, 1 * sizeof(cuDoubleComplex), cudaMemcpyDeviceToHost);

    //printf("GPU result : %18.12f %18.12f\n", h_zresult[0].x, h_zresult[0].y);

    //cudaFree(d_zresult);
    //free(h_zresult);
    //free(zs2);
}

template void spline_inner_product_gpu<double>(int size, double* r_dr, double* s1_coefs, double* s2_coefs);
