     gsl_matrix *Q = gsl_matrix_alloc(p_basis->num, p_basis->num);
     form_Q(p_basis, Q);
     */
+
+    t1 = clock();
+    time_in_sec = (t1 - t0) / (double)CLOCKS_PER_SEC;
+    time_txt += "Time_SCF_Init = " + std::to_string(time_in_sec) + " sec\n";
+    time_total += time_in_sec;
 
 
     // start SCF iterations


int main(int argc, char* argv[])
{
    // use spherical harmonic d function
    const int use_5d = 1;


         dev_xb,dev_yb,dev_zb, dev_lb,dev_mb,dev_nb, dev_bexps,dev_bcoef,
         n_combi, dev_start_contr, dev_end_contr, dev_eri);

    t1 = clock();
    time_in_sec = (t1 - t0) / (double)CLOCKS_PER_SEC;
    time_txt += "Time_2e_GPU   = " + std::to_string(time_in_sec) + " sec\n";
    time_total += time_in_sec;


    // copy the results back to host
    my_cuda_safe(cudaMemcpy(h_eri, dev_eri, n_ERI_bytes, cudaMemcpyDeviceToHost),"mem_eri"); 

    t0 = clock();
    time_in_sec = (t0 - t1) / (double)CLOCKS_PER_SEC;
    time_txt += "Time_2e_HMem  = " + std::to_string(time_in_sec) + " sec\n";
    time_total += time_in_sec;


    /* just for test...
