diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 14be0eb..82b0f85 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -7,32 +7,54 @@ target_include_directories(sgpu.exe PRIVATE
     "${CMAKE_SOURCE_DIR}/src/libcintw"
     "${CMAKE_SOURCE_DIR}/include"
     )
+
 set(TEST_COMP_FLAGS ${CMAKE_CXX_FLAGS})
+
+# CUDA/GPU dependencies - only when USE_ACC is enabled
 if(USE_ACC)
+    enable_language(CUDA)
+    find_package(CUDA REQUIRED)
+    find_library(CUDART_LIBRARY cudart ${CMAKE_CUDA_IMPLICIT_LINK_DIRECTORIES})
+    
     target_link_options(sgpu.exe PRIVATE 
         -acc=gpu
         )
     set(TEST_COMP_FLAGS ${TEST_COMP_FLAGS} 
         -DUSE_ACC=1 -acc=gpu
         )
+    
+    # GPU-specific link libraries
+    set(GPU_LIBRARIES 
+        "${CUDART_LIBRARY}" 
+        "${CUDA_LIBRARIES}" 
+        "${CUDA_TOOLKIT_ROOT_DIR}/../../math_libs/lib64/libcublas.so" 
+        "${CUDA_TOOLKIT_ROOT_DIR}/../../math_libs/lib64/libcusolver.so"
+        )
+else()
+    # CPU-only mode
+    set(GPU_LIBRARIES "")
+    message(STATUS "Building CPU-only version (USE_ACC=OFF)")
 endif()
 
 #LAPACK cmake not right.
 #need to use this as a variable below
 find_package(LAPACK)
-
 find_package(BLAS)
-enable_language(CUDA)
-find_package(CUDA REQUIRED)
-find_library(CUDART_LIBRARY cudart ${CMAKE_CUDA_IMPLICIT_LINK_DIRECTORIES})
-
-#enable_language(FORTRAN)
-#find_library(FORTRANLIBS fortranlibs ${CMAKE_CUDA_IMPLICIT_LINK_DIRECTORIES})
 
 set_property(TARGET sgpu.exe PROPERTY CXX_STANDARD 14)
 target_compile_options(sgpu.exe PRIVATE ${TEST_COMP_FLAGS} -O2)
-#target_link_libraries(sgpu.exe PRIVATE SlaterGPU io "${CUDART_LIBRARY}" "${CUDA_LIBRARIES}" )
-target_link_libraries(sgpu.exe PRIVATE SlaterGPU io cintw "${FORTRANLIBS}" "${CUDART_LIBRARY}" "${CUDA_LIBRARIES}" "${CUDA_TOOLKIT_ROOT_DIR}/../../math_libs/lib64/libcublas.so" "${CUDA_TOOLKIT_ROOT_DIR}/../../math_libs/lib64/libcusolver.so" "-L/home/paulzim/integrate/gpu/lib/lapack-3.9.0/lib/ -llapack -lblas" "-L/export/apps/RockyOS8/nvidia_hpc/2024_249/Linux_x86_64/24.9/ -fortranlibs" )
+
+# Link libraries - GPU_LIBRARIES will be empty for CPU-only builds
+target_link_libraries(sgpu.exe PRIVATE 
+    SlaterGPU 
+    io 
+    cintw 
+    "${FORTRANLIBS}" 
+    ${GPU_LIBRARIES}
+    "-L/home/paulzim/integrate/gpu/lib/lapack-3.9.0/lib/ -llapack -lblas" 
+    "-L/export/apps/RockyOS8/nvidia_hpc/2024_249/Linux_x86_64/24.9/ -fortranlibs"
+    )
+
 file(COPY lih_VK1 DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 file(COPY lih_ri5 DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 file(COPY h2_ri5 DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
@@ -42,6 +64,5 @@ file(COPY noplus_ri5 DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 file(COPY ch4_ri5 DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 file(COPY run.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 
-# Add installation for the executable
 install(TARGETS sgpu.exe
-        RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})
+        RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})
\ No newline at end of file
diff --git a/examples/main.cpp b/examples/main.cpp
index 26baeab..d48726d 100644
--- a/examples/main.cpp
+++ b/examples/main.cpp
@@ -68,14 +68,20 @@ void compute_ps_integrals_to_disk(int natoms, int* atno, double* coords, vector<
   for (int j=0;j<N2;j++) Ssp[j] = Tsp[j] = Ensp[j] = pVpsp[j] = 0.;
 
   compute_STEn_ps(natoms,atno,coords,basis,nquad,nmu,nnu,nphi,Ssp,Tsp,Ensp,prl);
-  //compute_pVp_ps(natoms,atno,coords,basis,nquad,nmu,nnu,nphi,pVpsp,prl);
-  //compute_pVp_3c_ps(natoms,atno,coords,basis,nquad,nquad2,nsplit,nmu,nnu,nphi,pVpsp,prl);
+ 
+  bool x2c = read_int("X2C");
+  if (x2c)
+  {
+    compute_pVp_ps(natoms,atno,coords,basis,nquad,nmu,nnu,nphi,pVpsp,prl);
+    compute_pVp_3c_ps(natoms,atno,coords,basis,nquad,nquad2,nsplit,nmu,nnu,nphi,pVpsp,prl);
+  }
   compute_2c_ps(0,0,gamma,nbatch,natoms,atno,coords,basis_aux,nquad,nmu,nnu,nphi,Asp,prl);
   compute_3c_ps(0,0,gamma,nbatch,natoms,atno,coords,basis,basis_aux,nquad,nquad2,nsplit,nmu,nnu,nphi,Ensp,Csp,prl);
 
   if (prl > 0) printf("Printing PS Integral Files:\n");
   write_S_En_T(N,Ssp,Ensp,Tsp);
-  //write_square(N,pVpsp,"pVp",2);
+  if (x2c)
+    write_square(N,pVpsp,"pVp",2);
   write_square(Naux,Asp,"A",2);
   write_C(Naux,N2,Csp);
 
@@ -297,7 +303,9 @@ int main(int argc, char* argv[]) {
       if (prl > 0) printf("Printing Standard Integral Files:\n");
       write_S_En_T(N,S,En,T);
       write_square(Naux,A,"A",2);
-      write_square(N,pVp,"pVp",2);
+      bool x2c = read_int("X2C");
+      if (x2c)
+        write_square(N,pVp,"pVp",2);
       write_C(Naux, N2, C);
 
       if(c4 > 0)
diff --git a/include/cintwrapper.h b/include/cintwrapper.h
index 8e6d642..041f835 100644
--- a/include/cintwrapper.h
+++ b/include/cintwrapper.h
@@ -99,6 +99,10 @@ void get_overlap_ri(double * overlap, int Naux,
                     int natm, int nbas, int nbas_ri, int nenv,
                     int *atm, int* bas, double *env);
 
+void get_hcore(double *hcore, int N,
+               int natm, int nbas, int nenv,
+               int *atm, int *bas, double *env);
+
 void get_hcore(double *hcore, double *En, double *T, int N,
                int natm, int nbas, int nenv,
                int *atm, int *bas, double *env);
diff --git a/include/cpu_util.h b/include/cpu_util.h
index 9a259d7..0fe6f21 100644
--- a/include/cpu_util.h
+++ b/include/cpu_util.h
@@ -52,6 +52,6 @@ void mat_times_mat_bt_cpu(double* C, double* A, double* B, int M, int N, int K,
 void mat_times_mat_bt_cpu(float* C, float* A, float* B, int N);
 void mat_times_mat_bt_cpu(double* C, double* A, double* B, int M, int N, int K);
 void mat_times_mat_bt_cpu(double* C, double* A, double* B, int N);
-
+void mat_times_mat_at_cpu(double* C, double* A, double* B, int M, int N, int K);
 
 #endif
diff --git a/include/cuda_util.h b/include/cuda_util.h
index 7e35284..3f58030 100644
--- a/include/cuda_util.h
+++ b/include/cuda_util.h
@@ -1,9 +1,16 @@
 #ifndef CUDA_UTILH
 #define CUDA_UTILH
 
-//#include <cublas.h>
+#if !USE_ACC
+#include <complex>
+typedef std::complex<double> cuDoubleComplex;
+typedef std::complex<float> cuFloatComplex;
+typedef void* cusolverDnHandle_t;
+typedef void* cublasHandle_t;
+#else
 #include <cublas_v2.h>
 #include <cusolverDn.h>
+#endif
 
 int invert_eigen_cusolver(int size, double* A, double eig_max, cusolverDnHandle_t& cu_hdl);
 int invert_stable_cusolver(int size, double* A, double delta, cusolverDnHandle_t& cu_hdl);
@@ -41,4 +48,5 @@ void copy_to_all_gpu(int ngpu, int s1, double* A, int include_first);
 void copy_to_all_gpu(int ngpu, int s1, float* A, int include_first);
 void copy_to_all_gpu(int ngpu, int s1, int s2, double** A, int include_first);
 
-#endif
+
+#endif // CUDA_UTILH
\ No newline at end of file
diff --git a/include/cusp.h b/include/cusp.h
index 2b070b7..fcefc37 100644
--- a/include/cusp.h
+++ b/include/cusp.h
@@ -24,6 +24,7 @@
 
 using namespace std;
 
+
 void compute_diatomic_symm(int natoms, int* atno, vector<vector<double> > basis, vector<double*>& pB_all, int prl);
 
 void compute_cusp(int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* pB1, double* pB2, int prl);
diff --git a/include/integrals.h b/include/integrals.h
index b475b07..b6a0a57 100644
--- a/include/integrals.h
+++ b/include/integrals.h
@@ -228,6 +228,7 @@ void copy_symm_4c_ps_cpu(int natoms, int* n2i, int N, double* olp);
 int get_natoms_with_basis(int natoms, int* atno, vector<vector<double> >& basis);
 
 vector<vector<double> > setup_integrals_gsgpu(vector<vector<double> >& basis_aux, int natoms, int* atno, double* coords, int& nbas, int& nenv, int& N, int& Naux, int& nbas_ri, int* &atm, int* &bas, double* &env, int prl);
+void compute_integrals_g(int natm, int nbas, int nenv, int N, int Naux, int nbas_ri, int* atm, int* bas, double* env, double* S, double* T, double* jH1, double* A, double* C, int prl);
 void compute_integrals_g(int natm, int nbas, int nenv, int N, int Naux, int nbas_ri, int* atm, int* bas, double* env, double* S, double* En, double* T, double* jH1, double* A, double* C, double* pvp, int prl);
 void compute_gaussian_integrals_to_disk(int N, int Naux, int natoms, int nbas, int nenv, int nbas_ri, int* atm, int* bas, double* env);
 
diff --git a/src/integrals/CMakeLists.txt b/src/integrals/CMakeLists.txt
index c11e217..0c0fcc3 100644
--- a/src/integrals/CMakeLists.txt
+++ b/src/integrals/CMakeLists.txt
@@ -1,4 +1,5 @@
-set(INTEGRALS_SOURCES
+# Common sources for both CPU and GPU builds
+set(INTEGRALS_COMMON_SOURCES
     integrals.cpp
     integrals_d.cpp
     integrals_aux.cpp
@@ -19,17 +20,27 @@ set(INTEGRALS_SOURCES
     ps_first_second_order.cpp
     cusp.cpp
     cpu_util.cpp
-    gpu_util.cpp
     grid_util.cpp
-    cuda_util.cpp
     scf_util.cpp
     symm.cpp
     gauss.cpp
     hess.cpp
     opt.cpp
     jellium_ints.cpp
+    gpu_util.cpp
+    cuda_util.cpp
     )
 
+
+# Combine sources based on USE_ACC flag
+if(USE_ACC)
+    set(INTEGRALS_SOURCES ${INTEGRALS_COMMON_SOURCES})
+    message(STATUS "Building SlaterGPU library with GPU support")
+else()
+    set(INTEGRALS_SOURCES ${INTEGRALS_COMMON_SOURCES})
+    message(STATUS "Building SlaterGPU library for CPU-only")
+endif()
+
 add_library(SlaterGPU STATIC "${INTEGRALS_SOURCES}")
 target_include_directories(SlaterGPU PUBLIC
     $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/include>
@@ -52,4 +63,4 @@ else()
     target_compile_options(SlaterGPU PRIVATE
         ${SLATER_COMP_FLAGS} -DUSE_ACC=0
         )
-endif()
+endif()
\ No newline at end of file
diff --git a/src/integrals/cpu_util.cpp b/src/integrals/cpu_util.cpp
index 44b39dc..7afd923 100644
--- a/src/integrals/cpu_util.cpp
+++ b/src/integrals/cpu_util.cpp
@@ -1120,3 +1120,23 @@ void cross(double* m, double* r1, double* r2)
 
   return;
 }
+
+// Extended version of mat_times_mat_at_cpu with M, N, K parameters
+void mat_times_mat_at_cpu(double* C, double* A, double* B, int M, int N, int K)
+{
+  // C = A^T * B
+  // A is K x M, B is K x N, C is M x N
+  
+  int LDA = M;
+  int LDB = N;
+  int LDC = N;
+
+  double ALPHA = 1.0;
+  double BETA = 0.0;
+
+  char TA = 'T';
+  char TB = 'N';
+  dgemm_(&TB,&TA,&N,&M,&K,&ALPHA,B,&LDB,A,&LDA,&BETA,C,&LDC);
+
+  return;
+}
\ No newline at end of file
diff --git a/src/integrals/cuda_util.cpp b/src/integrals/cuda_util.cpp
index 98327e0..d3c6763 100644
--- a/src/integrals/cuda_util.cpp
+++ b/src/integrals/cuda_util.cpp
@@ -1,3 +1,5 @@
+//Vaibhav needs to fix this
+
 #include <stdio.h>
 #include <stdlib.h>
 #include <unistd.h>
@@ -10,7 +12,10 @@
 #include "cuda_util.h"
 #include "cpu_util.h"
 
+#if USE_ACC
 #include "cuda_runtime.h"
+#include <cuComplex.h>
+#endif
 
 //these functions assume the matrices are already on the device, via ACC
 //hoping to reduce overhead by getting cusolver handle just once
@@ -23,21 +28,32 @@ void print_square_sm(int N, float* A);
 void print_square_sm(int N, double* A);
 double read_float(string filename);
 
-#include <cuComplex.h>
+
 
 void print_square_complex(int N, cuDoubleComplex* A)
 {
+ #if !USE_ACC
+  printf("\n ERROR: there is no CPU implementation for print_square_complex \n"); exit(-1);
+ #else
   for (int i=0;i<N;i++)
   {
     for (int j=0;j<N;j++)
       printf("  %8.5f + %8.5fi",A[i*N+j].x,A[i*N+j].y);
     printf("\n");
   }
+
+ #endif
 }
 
 //void expmat_complex(int N, double* theta1, double* theta1i, double* U, cusolverDnHandle_t cu_hdl, cublasHandle_t cublasH)
 double expmat_complex(int N, double* theta1, double* theta1i, double* jCA, double* U, cusolverDnHandle_t cu_hdl, cublasHandle_t cublasH)
 {
+  double ss = 0.;
+ 
+  #if !USE_ACC
+    printf("\n ERROR: there is no CPU implementation for expmat_complex \n"); exit(-1);
+    return 0.;
+  #else
  //this is a bit of a mess
  //takes exp(-i theta)
 
@@ -140,8 +156,7 @@ double expmat_complex(int N, double* theta1, double* theta1i, double* jCA, doubl
   cudaDeviceSynchronize();
 
  //CPMZ cannot use this in current form
-  double ss = 0.;
-  for (int i=0;i<N;i++)
+ for (int i=0;i<N;i++)
   for (int j=0;j<i;j++)
   {
     double x1 = U1[i*N+j].x;
@@ -194,6 +209,7 @@ double expmat_complex(int N, double* theta1, double* theta1i, double* jCA, doubl
   delete [] Ue;
   delete [] eix;
   delete [] U1;
+ #endif
 
   return ss;
 }
@@ -203,6 +219,11 @@ void expmat(int N, double* theta1, double* U, cusolverDnHandle_t cu_hdl, cublasH
  //this is a bit of a mess
  //takes exp(-i theta)
 
+ #if !USE_ACC
+  printf("\n ERROR: there is no CPU implementation for expmat \n"); exit(-1);
+  return;
+ #else
+
   int N2 = N*N;
 
   cusolverEigMode_t jobz = CUSOLVER_EIG_MODE_VECTOR;
@@ -304,11 +325,16 @@ void expmat(int N, double* theta1, double* U, cusolverDnHandle_t cu_hdl, cublasH
   delete [] eix;
   delete [] U1;
 
+ #endif
+
   return;
 }
 
 void mat_times_mat(float* C, float* A, float* B, int M, int N, int K)
 {
+  #if !USE_ACC
+    return mat_times_mat_cpu(C,A,B,M,N,K);
+  #endif
   int MN = M*N;
   int MK = M*K;
   int NK = N*K;
@@ -336,6 +362,9 @@ void mat_times_mat(float* C, float* A, float* B, int M, int N, int K)
 
 void mat_times_mat(float* C, float* A, float* B, int N)
 {
+  #if !USE_ACC
+     return mat_times_mat_cpu(C,A,B,N);
+  #endif
   return mat_times_mat(C,A,B,N,N,N);
 }
 
@@ -612,6 +641,9 @@ int mat_root_inv_cusolver(double* A, int size, cusolverDnHandle_t& cu_hdl)
 
 void diagonalize_cusolver(int N, float* A, float* Ae, cusolverDnHandle_t& cu_hdl)
 {
+ #if !USE_ACC
+ printf("\n ERROR: there is no CPU implementation for float diagonalize_cusolver \n"); exit(-1);
+ #else
   cusolverEigMode_t jobz = CUSOLVER_EIG_MODE_VECTOR;
   cublasFillMode_t uplo = CUBLAS_FILL_MODE_LOWER;
 
@@ -632,7 +664,7 @@ void diagonalize_cusolver(int N, float* A, float* Ae, cusolverDnHandle_t& cu_hdl
 
  #pragma acc exit data delete(work[0:Lwork],info[0:1])
 
-#if 0
+ #if 0
   #pragma acc update self(A[0:N*N],Ae[0:N])
   printf("  eigenvalues: ");
   for (int i=0;i<N;i++)
@@ -640,16 +672,21 @@ void diagonalize_cusolver(int N, float* A, float* Ae, cusolverDnHandle_t& cu_hdl
   printf("\n");
   printf("  eigenvectors: \n");
   print_square_sm(N,A);
-#endif
+ #endif
 
   delete [] work;
   delete [] info;
 
+ #endif
+
   return;
 }
 
 void diagonalize_cusolver(int Ne, int N, double* A, double* Ae, cusolverDnHandle_t& cu_hdl)
 {
+ #if !USE_ACC
+  return la_diag(N,N,A,Ae);
+ #else
   cusolverEigMode_t jobz = CUSOLVER_EIG_MODE_VECTOR;
   cublasFillMode_t uplo = CUBLAS_FILL_MODE_LOWER;
 
@@ -699,6 +736,7 @@ void diagonalize_cusolver(int Ne, int N, double* A, double* Ae, cusolverDnHandle
   delete [] work;
   delete [] info;
 
+ #endif
   return;
 }
 
@@ -936,6 +974,11 @@ int invert_cusolver(int N, float* A, cusolverDnHandle_t& cu_hdl)
 {
   int reset = 0;
 
+ #if !USE_ACC
+  printf("\n ERROR: there is no CPU implementation for invert_cusolver \n"); exit(-1);
+  return reset;
+ #else
+
   int N2 = N*N;
   float* B = new float[N2];
   #pragma acc enter data create(B[0:N2])
@@ -1015,12 +1058,19 @@ int invert_cusolver(int N, float* A, cusolverDnHandle_t& cu_hdl)
   delete [] info;
   delete [] work;
 
+ #endif
+
   return reset;
 }
 
 int invert_cusolver(int N, double* A, cusolverDnHandle_t& cu_hdl)
 {
   int reset = 0;
+ #if !USE_ACC
+  printf("\n ERROR: there is no CPU implementation for invert_cusolver \n"); exit(-1);
+  return reset;
+ #else
+
 
   int N2 = N*N;
   double* B = new double[N2];
@@ -1114,6 +1164,7 @@ int invert_cusolver(int N, double* A, cusolverDnHandle_t& cu_hdl)
   cudaFree(info);
   cudaFree(work);
 
+ #endif
   return reset;
 }
 
diff --git a/src/integrals/gauss.cpp b/src/integrals/gauss.cpp
index 7be6b93..7bbdc27 100644
--- a/src/integrals/gauss.cpp
+++ b/src/integrals/gauss.cpp
@@ -240,11 +240,13 @@ void integrate_hole_para_gh(double* rdm, bool full_rdm, bool hfx_on, int Nc, int
     val2b[i1][j] = v2;
   }
 
+  #if USE_ACC
   copy_to_all_gpu(ngpu,gsa6,grid,0);
   for (int i1=0;i1<N;i1++)
     copy_to_all_gpu(ngpu,gsa,val2[i1],0);
   for (int i1=0;i1<N;i1++)
     copy_to_all_gpu(ngpu,gsba,val2b[i1],0);
+  #endif
 
 
   int gdiv = gsba/ngpu+1;
diff --git a/src/integrals/gpu_util.cpp b/src/integrals/gpu_util.cpp
index 2e2c40c..965780a 100644
--- a/src/integrals/gpu_util.cpp
+++ b/src/integrals/gpu_util.cpp
@@ -3,6 +3,8 @@
 
 void copy_to_all_gpu(int ngpu, int s1, double* A, int include_first)
 {
+ #if USE_ACC
+
   if (ngpu<1) return;
   if (include_first==0)
   {
@@ -21,11 +23,16 @@ void copy_to_all_gpu(int ngpu, int s1, double* A, int include_first)
   }
 
   acc_set_device_num(0,acc_device_nvidia);
+
+ #endif
+
   return;
 }
 
 void copy_to_all_gpu(int ngpu, int s1, float* A, int include_first)
 {
+ #if USE_ACC
+
   if (ngpu<1) return;
   if (include_first==0)
   {
@@ -44,11 +51,16 @@ void copy_to_all_gpu(int ngpu, int s1, float* A, int include_first)
   }
 
   acc_set_device_num(0,acc_device_nvidia);
+
+ #endif
+
   return;
 }
 
 void copy_to_all_gpu(int ngpu, int s1, int s2, double** A, int include_first)
 {
+ #if USE_ACC
+
   if (ngpu<1) return;
   if (include_first==0)
   {
@@ -67,5 +79,8 @@ void copy_to_all_gpu(int ngpu, int s1, int s2, double** A, int include_first)
   }
 
   acc_set_device_num(0,acc_device_nvidia);
+
+ #endif
+  
   return;
 }
\ No newline at end of file
diff --git a/src/integrals/integrals_aux.cpp b/src/integrals/integrals_aux.cpp
index caf9cea..4a52e24 100644
--- a/src/integrals/integrals_aux.cpp
+++ b/src/integrals/integrals_aux.cpp
@@ -5,6 +5,7 @@
 #include "../../include/cintprep.h"
 
 #include "write.h"
+#include "read.h"
 
 void auto_crash()
 {
@@ -1617,13 +1618,64 @@ vector<vector<double> > setup_integrals_gsgpu(vector<vector<double> >& basis_aux
   return basis;
 }
 
+void compute_integrals_g(int natm, int nbas, int nenv, int N, int Naux, int nbas_ri, int* atm, int* bas, double* env, double* S, double* T, double* jH1, double* A, double* C, int prl)
+{
+  get_overlap(S, N, natm, nbas, nenv, atm, bas, env);
+  get_hcore(jH1, N, natm, nbas, nenv, atm, bas, env);
+  if (T!=NULL)
+    get_tcore(T, N, natm, nbas, nenv, atm, bas, env);
+
+  if (prl>1)
+  {
+    printf("\n S: \n");
+    print_square(N,S);
+    printf("\n jH1: \n");
+    print_square(N,jH1);
+  }
+
+  if (Naux>0)
+  {
+    gen_eri_2c(A, Naux, natm, nbas, nenv, nbas_ri, atm, bas, env);
+    gen_eri_3c(C, N, Naux, natm, nbas, nenv, nbas_ri, atm, bas, env);
+  }
+  else
+  {
+    printf("  no auxiliary basis \n");
+  }
+
+  if (prl>1)
+  {
+    printf("\n A: \n");
+    for (int m=0;m<Naux;m++)
+    {
+      for (int n=0;n<Naux;n++)
+        printf(" %8.5f",A[m*Naux+n]);
+      printf("\n");
+    }
+
+    int N2 = N*N;
+    printf("\n C: \n");
+    for (int m=0;m<N2;m++)
+    {
+      for (int n=0;n<Naux;n++)
+        printf(" %8.5f",C[m*Naux+n]);
+      printf("\n");
+    }
+  }
+
+  return;
+}
+
 void compute_integrals_g(int natm, int nbas, int nenv, int N, int Naux, int nbas_ri, int* atm, int* bas, double* env, double* S, double* En, double* T, double* jH1, double* A, double* C, double* pvp, int prl)
 {
   get_overlap(S, N, natm, nbas, nenv, atm, bas, env);
   get_hcore(jH1, En, T, N, natm, nbas, nenv, atm, bas, env);
   //if (T!=NULL)
   //  get_tcore(T, N, natm, nbas, nenv, atm, bas, env);
-  gen_pvp(pvp, N, natm, nbas, nenv, atm, bas, env);
+
+  bool x2c = read_int("X2C");
+  if (x2c)
+    gen_pvp(pvp, N, natm, nbas, nenv, atm, bas, env);
 
   if (prl>1)
   {
diff --git a/src/integrals/integrals_d.cpp b/src/integrals/integrals_d.cpp
index c9cc065..c494e99 100644
--- a/src/integrals/integrals_d.cpp
+++ b/src/integrals/integrals_d.cpp
@@ -8,7 +8,9 @@
 #include <string>
 void auto_crash();
 void print_duration(chrono::high_resolution_clock::time_point t1, chrono::high_resolution_clock::time_point t2, string name);
-void copy_to_all_gpu(int ngpu, int s1, double* A, int include_first);
+#if USE_ACC
+  void copy_to_all_gpu(int ngpu, int s1, double* A, int include_first);
+#endif
 
 void gather_12_d_En_0(int s1, int s2, int gs, int iN, float** valtx1, float** valtx2, float** valS1x, float** valS2x, float* wtt1, float* wtt2)
 {
@@ -292,7 +294,9 @@ void compute_d_3c_para(int npgu, int natoms, int* atno, float* coords, vector<ve
   }
 
  //distribute dC to all gpus
+ #if USE_ACC
   copy_to_all_gpu(nomp,N2a,dC,0);
+ #endif
 
 #if USE_ACC
  #pragma omp parallel for num_threads(nomp)
diff --git a/src/integrals/prosph.cpp b/src/integrals/prosph.cpp
index f064e6a..4e8ff02 100644
--- a/src/integrals/prosph.cpp
+++ b/src/integrals/prosph.cpp
@@ -10,8 +10,9 @@
 // 2. //asynchronous?
 
 // get_cfn turned off
-
-void copy_to_all_gpu(int ngpu, int s1, double* A, int include_first);
+#if USE_ACC
+  void copy_to_all_gpu(int ngpu, int s1, double* A, int include_first);
+#endif
 void swap_grids_gpu(int gs, double*& grid, double*& gridm, double*& wt, double* gridn, double* gridmn, double* wtn);
 void auto_crash();
 
@@ -1501,7 +1502,9 @@ void do_3c_integrals_ps(const double epsilon, double cf, int nomp, int ss, int m
     int nphi1 = nphi;
 
     get_3c_position(A1,B1,C1,A2,B2,C2,A3,B3,C3,rot);
-    copy_to_all_gpu(nomp,9,rot,1);
+    #if USE_ACC
+      copy_to_all_gpu(nomp,9,rot,1);
+    #endif
     double coordn[9];
     coordn[0] = A1; coordn[1] = B1; coordn[2] = C1;
     coordn[3] = A2; coordn[4] = B2; coordn[5] = C2;
@@ -1708,7 +1711,9 @@ void do_3c_integrals_ps(const double epsilon, double cf, int nomp, int ss, int m
     int nphi1 = nphi;
 
     get_3c_position(A1,B1,C1,A2,B2,C2,A3,B3,C3,rot);
-    copy_to_all_gpu(nomp,9,rot,1);
+    #if USE_ACC
+      copy_to_all_gpu(nomp,9,rot,1);
+    #endif
     double coordn[9];
     coordn[0] = A1; coordn[1] = B1; coordn[2] = C1;
     coordn[3] = A2; coordn[4] = B2; coordn[5] = C2;
@@ -2204,7 +2209,9 @@ void integrate_STEnAC_2c(int natoms, int* atno, double* coords, vector<vector<do
     {
       do_rot = 1;
       get_2c_position(A1,B1,C1,A2,B2,C2,rot);
-      copy_to_all_gpu(nomp,9,rot,1);
+      #if USE_ACC
+        copy_to_all_gpu(nomp,9,rot,1);
+      #endif
     }
     else
       A1 = B1 = C1 = A2 = B2 = C2 = 0.;
diff --git a/src/integrals/reduce.cpp b/src/integrals/reduce.cpp
index 01efb8a..9b297b8 100644
--- a/src/integrals/reduce.cpp
+++ b/src/integrals/reduce.cpp
@@ -354,9 +354,9 @@ void reduce_2c1(int tid, int s1, int s2, int gs, double** val1, double** val3, i
 
 #else
   for (int i1=s1;i1<s2;i1++)
-  for (int i2=s3;i2<s4;i2++)
+  for (int i2=s1;i2<s2;i2++)
   {
-    int ii1 = i1-s1; int ii2 = i2-s3;
+    int ii1 = i1-s1; int ii2 = i2-s1;
     double val = 0.;
     for (int j=0;j<gs;j++)
       val += val1[ii1][j]*val3[ii2][j];
diff --git a/src/integrals/scf_util.cpp b/src/integrals/scf_util.cpp
index 69a1990..63dc3b2 100644
--- a/src/integrals/scf_util.cpp
+++ b/src/integrals/scf_util.cpp
@@ -223,10 +223,10 @@ void ao_to_mo(int N, double* jAO, double* jMO, double* jCA, double* tmp)
 {
  #if !USE_ACC
   return ao_to_mo_cpu(N,jAO,jMO,jCA,tmp);
- #endif
+ #else
   mat_times_mat(tmp,jAO,jCA,N);
   mat_times_mat_at(jMO,jCA,tmp,N);
-
+#endif
   return;
 }
 
diff --git a/src/libcintw/cintwrapper.cpp b/src/libcintw/cintwrapper.cpp
index 444fd9e..d9a82ea 100644
--- a/src/libcintw/cintwrapper.cpp
+++ b/src/libcintw/cintwrapper.cpp
@@ -157,6 +157,127 @@ void get_overlap_ri(double * overlap, int Naux,
   } // else
 }
 
+void get_hcore(double *hcore, int N,
+               int natm, int nbas, int nenv,
+               int *atm, int *bas, double *env) {
+  int idx_i = 0;
+  int di, dj;
+  int shls[2];
+
+  double *tmp = new double[N * N];
+  
+  if (BT::DO_CART) {
+    for (int i = 0; i < nbas; i++) {
+      int idx_j = 0;
+      shls[0] = i;
+      di = CINTcgto_cart(i, bas);
+      for (int j = 0; j <= i; j++) {
+        dj = CINTcgto_cart(j, bas);
+        shls[1] = j;
+        double *buf = new double[di*dj];
+        cint1e_nuc_cart(buf,shls,atm,natm,bas,nbas,env);
+
+        for (int j1 = 0; j1 < dj; j1++) {
+          for (int i1 = 0; i1 < di; i1++) {
+            int oi = i1 + idx_i;
+            int oj = j1 + idx_j;
+            tmp[oi * N + oj] = tmp[oj * N + oi] = buf[j1*di + i1];
+          } // for i1
+        } // for j1
+        delete [] buf;
+        idx_j += dj;
+      } // for j
+      idx_i += di;
+    }// for i
+
+    for (int i = 0; i < N * N; i++) {
+      hcore[i] = tmp[i];
+    }
+
+    idx_i = 0;
+    for (int i = 0; i < nbas; i++) {
+      int idx_j = 0;
+      shls[0] = i;
+      di = CINTcgto_cart(i, bas);
+      for (int j = 0; j <= i; j++) {
+        dj = CINTcgto_cart(j, bas);
+        shls[1] = j;
+        double *buf = new double[di*dj];
+        cint1e_kin_cart(buf,shls,atm,natm,bas,nbas,env);
+
+        for (int j1 = 0; j1 < dj; j1++) {
+          for (int i1 = 0; i1 < di; i1++) {
+            int oi = i1 + idx_i;
+            int oj = j1 + idx_j;
+            tmp[oi * N + oj] = tmp[oj * N + oi] = buf[j1*di + i1];
+          } // for i1
+        } // for j1
+        delete [] buf;
+        idx_j += dj;
+      } // for j
+      idx_i += di;
+    }// for i
+  } // if BT::DO_CART
+  else {
+    for (int i = 0; i < nbas; i++) {
+      int idx_j = 0;
+      shls[0] = i;
+      di = CINTcgto_spheric(i, bas);
+      for (int j = 0; j <= i; j++) {
+        dj = CINTcgto_spheric(j, bas);
+        shls[1] = j;
+        double *buf = new double[di*dj];
+        cint1e_nuc_sph(buf,shls,atm,natm,bas,nbas,env);
+
+        for (int j1 = 0; j1 < dj; j1++) {
+          for (int i1 = 0; i1 < di; i1++) {
+            int oi = i1 + idx_i;
+            int oj = j1 + idx_j;
+            tmp[oi * N + oj] = tmp[oj * N + oi] = buf[j1*di + i1];
+          } // for i1
+        } // for j1
+        delete [] buf;
+        idx_j += dj;
+      } // for j
+      idx_i += di;
+    }// for i
+
+    for (int i = 0; i < N * N; i++) {
+      hcore[i] = tmp[i];
+    }
+
+    idx_i = 0;
+    for (int i = 0; i < nbas; i++) {
+      int idx_j = 0;
+      shls[0] = i;
+      di = CINTcgto_spheric(i, bas);
+      for (int j = 0; j <= i; j++) {
+        dj = CINTcgto_spheric(j, bas);
+        shls[1] = j;
+        double *buf = new double[di*dj];
+        cint1e_kin_sph(buf,shls,atm,natm,bas,nbas,env);
+
+        for (int j1 = 0; j1 < dj; j1++) {
+          for (int i1 = 0; i1 < di; i1++) {
+            int oi = i1 + idx_i;
+            int oj = j1 + idx_j;
+            tmp[oi * N + oj] = tmp[oj * N + oi] = buf[j1*di + i1];
+          } // for i1
+        } // for j1
+        delete [] buf;
+        idx_j += dj;
+      } // for j
+      idx_i += di;
+    }// for i
+  }  
+
+  for (int i = 0; i < N * N; i++) {
+    hcore[i] += tmp[i];
+  }
+
+  delete [] tmp;
+}
+
 void get_hcore(double *hcore, double *En, double *T, int N,
                int natm, int nbas, int nenv,
                int *atm, int *bas, double *env) {