diff --git a/include/becke.h b/include/becke.h
index 71b2362..fbfd156 100644
--- a/include/becke.h
+++ b/include/becke.h
@@ -27,7 +27,7 @@ void becke_weight_3c(int gs, float* grid1, float* wt1, float* grid2, float* wt2,
 void becke_weight_2c(int gs, float* grid1, float* wt1, float* grid2, float* wt2,
                      int Z1, int Z2, float A2, float B2, float C2);
 void becke_weight_2d(int gs, double* grid1, double* wt1, double* grid2, double* wt2,
-                     double zeta1, double zeta2, double A2, double B2, double C2);
+                     bool is_Z, double zeta1, double zeta2, double A2, double B2, double C2);
 
 float becke_ar(float r1, float r2);
 double becke_ard(double r1, double r2);
diff --git a/include/cint.h b/include/cint.h
index 57c78f2..ff013be 100644
--- a/include/cint.h
+++ b/include/cint.h
@@ -4,7 +4,8 @@
  * Parameters and function signature for libcint.
  */
 
-#define CINT_VERSION    4.4.0
+#define CINT_VERSION            "6.1.3"
+#define CINT_SOVERSION          6
 
 /* #undef I8 */
 #ifdef I8
@@ -49,9 +50,8 @@
 #define PTR_COORD       1
 #define NUC_MOD_OF      2
 #define PTR_ZETA        3
-#define PTR_FRAC_CHARGE 3
-#define RESERVE_ATMLOT1 4
-#define RESERVE_ATMLOT2 5
+#define PTR_FRAC_CHARGE 4
+#define RESERVE_ATMSLOT 5
 #define ATM_SLOTS       6
 
 
@@ -150,6 +150,78 @@ typedef struct {
     double **log_max_coeff;
     PairData **pairdata;  // NULL indicates not-initialized, NO_VALUE can be skipped
 } CINTOpt;
+
+// Add this macro def to make pyscf compatible with both v4 and v5
+#define HAVE_DEFINED_CINTENVVARS_H
+typedef struct {
+        FINT *atm;
+        FINT *bas;
+        double *env;
+        FINT *shls;
+        FINT natm;
+        FINT nbas;
+
+        FINT i_l;
+        FINT j_l;
+        FINT k_l;
+        FINT l_l;
+        FINT nfi;  // number of cartesian components
+        FINT nfj;
+        // in int1e_grids, the grids_offset and the number of grids
+        union {FINT nfk; FINT grids_offset;};
+        union {FINT nfl; FINT ngrids;};
+        FINT nf;  // = nfi*nfj*nfk*nfl;
+        FINT rys_order; // = nrys_roots for regular ERIs. can be nrys_roots/2 for SR ERIs
+        FINT x_ctr[4];
+
+        FINT gbits;
+        FINT ncomp_e1; // = 1 if spin free, = 4 when spin included, it
+        FINT ncomp_e2; // corresponds to POSX,POSY,POSZ,POS1, see cint.h
+        FINT ncomp_tensor; // e.g. = 3 for gradients
+
+        /* values may diff based on the g0_2d4d algorithm */
+        FINT li_ceil; // power of x, == i_l if nabla is involved, otherwise == i_l
+        FINT lj_ceil;
+        FINT lk_ceil;
+        FINT ll_ceil;
+        FINT g_stride_i; // nrys_roots * shift of (i++,k,l,j)
+        FINT g_stride_k; // nrys_roots * shift of (i,k++,l,j)
+        FINT g_stride_l; // nrys_roots * shift of (i,k,l++,j)
+        FINT g_stride_j; // nrys_roots * shift of (i,k,l,j++)
+        FINT nrys_roots;
+        FINT g_size;  // ref to cint2e.c g = malloc(sizeof(double)*g_size)
+
+        FINT g2d_ijmax;
+        FINT g2d_klmax;
+        double common_factor;
+        double expcutoff;
+        double rirj[3]; // diff by sign in different g0_2d4d algorithm
+        double rkrl[3];
+        double *rx_in_rijrx;
+        double *rx_in_rklrx;
+
+        double *ri;
+        double *rj;
+        double *rk;
+        // in int2e or int3c2e, the coordinates of the fourth shell
+        // in int1e_grids, the pointer for the grids coordinates
+        union {double *rl; double *grids;};
+
+        FINT (*f_g0_2e)();
+        void (*f_g0_2d4d)();
+        void (*f_gout)();
+        CINTOpt *opt;
+
+        /* values are assigned during calculation */
+        int *idx;
+        double ai[1];
+        double aj[1];
+        double ak[1];
+        double al[1];
+        double fac[1];
+        double rij[3];
+        double rkl[3];
+} CINTEnvVars;
 #endif
 
 FINT CINTlen_cart(const FINT l);
diff --git a/include/cintwrapper.h b/include/cintwrapper.h
index 91baf23..db0ffc9 100644
--- a/include/cintwrapper.h
+++ b/include/cintwrapper.h
@@ -141,9 +141,25 @@ void gen_4c_overlap(double* ovlp4, size_t N,
                     int natm, int nbas, int nenv,
                     int* atm, int* bas, double* env);
 
-void gen_eri(double **eri, int N,
+void gen_eri_4b(double** eri, int N1, int N2, int N3, int N,
+             int natm, int nbas1, int nbas2, int nbas3, int nbas, int nenv,
+             int* atm, int* bas, double* env, int nomp);
+
+void gen_eri_2b(double** eri, int N1, int N,
+             int natm, int nbas1, int nbas, int nenv,
+             int* atm, int* bas, double* env, int nomp);
+
+void gen_eri(double** eri, int N,
+             int natm, int nbas, int nenv,
+             int* atm, int* bas, double* env);
+
+void gen_eri(double** eri, int N,
+             int natm, int nbas, int nenv,
+             int* atm, int* bas, double* env, int nomp);
+
+void gen_eri_2(double** eri, int N,
              int natm, int nbas, int nenv,
-             int *atm, int *bas, double *env);
+             int* atm, int* bas, double* env, int nomp);
 
 void gen_jMOI_gto(double **eri, int N,
                   int natm, int nbas, int nenv,
diff --git a/include/cpu_util.h b/include/cpu_util.h
index a285d12..b9afbac 100644
--- a/include/cpu_util.h
+++ b/include/cpu_util.h
@@ -23,11 +23,15 @@ using namespace std;
 double randomf(double a, double b);
 
 void solve_axeb(int dim, double* A, double* b);
+void solve_axeb_ls(int m, int n, double* A, double* b);
+
+double determinant(int N, double* A);
 
 void expmat_complex_cpu(int N, double* theta, double* thetai, double* etheta);
 void expmat_cpu(int N, double *theta, double *etheta);
 
 int la_diagR(int neig, double* A, double* eigen, double* eigeni);
+void la_diag_gen(int neig, int s1, double* A, double* S, double* Ae);
 void la_diag(int neig, int s1, double* A, double* Ae);
 int invert_stable_cpu(double* A, int size, double delta);
 int invert_stable_cpu(double* A, int size, double delta, int prl);
diff --git a/include/cuda_util.h b/include/cuda_util.h
index 26b9b5b..c099cc6 100644
--- a/include/cuda_util.h
+++ b/include/cuda_util.h
@@ -22,6 +22,7 @@ void diagonalize_cusolver(int Ne, int N, double* A, double* Ae, cusolverDnHandle
 
 int mat_root_inv_cusolver(double* A, int size, cusolverDnHandle_t& cu_hdl);
 int mat_root_inv_stable_cusolver(double* A, int size, double delta, cusolverDnHandle_t& cu_hdl);
+int mat_root_inv_stable_cusolver(double* A, int size, double delta, cusolverDnHandle_t& cu_hdl, int prl);
 
 void solve_axeb_gpu(int dim, double* A, double* b, cusolverDnHandle_t& cusolverH);
 
@@ -51,4 +52,4 @@ 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 // CUDA_UTILH
\ No newline at end of file
+#endif // CUDA_UTILH
diff --git a/include/gauss.h b/include/gauss.h
index 571f9ec..3583d00 100644
--- a/include/gauss.h
+++ b/include/gauss.h
@@ -21,7 +21,7 @@ void eval_p_gh(int gs, float* grid, float* val, int n1, int l1, int m1, float no
 void eval_pd_gh(int gs, double* grid, double* val, int n1, int l1, int m1, double norm1, double zeta1);
 
 //Gaussian single-center potentials
-void eval_vghd(int gsa, double* grid, double* val, int n, int l, int m, double norm, double zt);
+void eval_vghd(int gsa, double* grid, double* val, int n, int l, double norm, double zt);
 
 void eval_hess_ghd(int gs, double* grid, double* val, int n1, int l1, int m1, double norm1, double zeta1);
 int eval_gh_full(int gs, float* grid, float** val1, int i1, int natoms, int nbas, int nenv, int N, int* atm, int* bas, double* env);
@@ -35,5 +35,7 @@ void wf_to_grid_gh(bool divide_by_rho, bool calc_rho, int natoms, int* atno, dou
 //void integrate_hole_para_gh(double* rdm, bool full_rdm, bool hfx_on, int Nc, int No, int M, int natoms, int* atno, double* coords, int gs, int gsb, vector<vector<double> >& basis,
 //                    double* Pao, double* Pmo, double* jCA, float* grid, float* gridb, float* wt, float* wtb, float* rho, float* vxch, int prl);
 
+void gen_4cj_overlap(double* ol4j, int nrad, int nang, double* ang_g, double* ang_w,
+                     int natoms, int* atno, double* coords, vector<vector<double> > basis, int prl);
 
 #endif
diff --git a/include/isotopes.h b/include/isotopes.h
new file mode 100644
index 0000000..5aaa9e6
--- /dev/null
+++ b/include/isotopes.h
@@ -0,0 +1,10 @@
+#ifndef ISOTOPESH
+#define ISOTOPESH
+
+//most abundant isotopes of elements H through Ne
+const int isos1[11] = { 0, 1, 4, 7, 9, 11, 12, 14, 16, 19, 20 };
+const int isos2[11] = { 0, 2, 3, 6, 0, 10, 13, 15, 18, 0, 22 };
+const int isos3[11] = { 0, 3, 0, 0, 0, 0, 0, 0, 17, 0, 21 };
+
+
+#endif
diff --git a/include/read.h b/include/read.h
index 4082568..e79e2c1 100644
--- a/include/read.h
+++ b/include/read.h
@@ -26,6 +26,7 @@ string SSTRF(float x);
 string SSTRF2(double x);
 
 vector<string> split1(const string &s, char delim);
+vector<string> split2(const string &s, char delim);
 
 int check_file(string filename);
 
@@ -33,6 +34,8 @@ int initialize(bool gbasis, vector<vector<double> >& basis, vector<vector<double
 bool get_secondary_basis(string name, int natoms, int* atno, double* coords, vector<vector<double> >& basis, vector<vector<double> >& basis_aux, int prl);
 void create_basis_aux(int natoms, vector<vector<double> >& basis_std, vector<vector<double> >& basis_aux);
 string read_basis_text(string aname);
+string read_basis_text(string filename, string aname);
+int read_geoms(int& ngeom1, int natoms, vector<double*>& coords, int prl);
 void print_basis(int natoms, vector<vector<double> >& basis, vector<vector<double> >& basis_aux, int prl);
 int read_rotate(int N, double* jCA);
 double read_float(string filename);
@@ -93,6 +96,10 @@ bool read_MOI_from_file(int Mm, int M3, double** MOI, string filename, int prl);
 
 bool check_PS();
 
+double nuclear_repulsion(int natoms, double* atno, double* coords);
+double nuclear_repulsion(int natoms, int* atno, double* coords);
+
+
 string get_aname(int Z);
 
 #endif
diff --git a/include/write.h b/include/write.h
index 6cee7a5..5856db0 100644
--- a/include/write.h
+++ b/include/write.h
@@ -51,6 +51,7 @@ void write_square_clean(int N, double* A, string fname, double thresh, int prl);
 
 void write_molden(bool gbasis, int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, double* eig, string fname);
 void write_molden(bool gbasis, int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, string fname);
+void write_molden_g(int btype, int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, double* eig, string fname);
 void write_molden_g(int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, double* eig, string fname);
 //void write_molden_vcf(int natoms, int* atno, double* coords, vector<vector<vcf> > &vcfs, string fname);
 void write_molden_ss(int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, string fname);
diff --git a/src/integrals/becke.cpp b/src/integrals/becke.cpp
index d489847..2305ca0 100644
--- a/src/integrals/becke.cpp
+++ b/src/integrals/becke.cpp
@@ -954,7 +954,7 @@ void becke_weight_2c(int gs, float* grid1, float* wt1, float* grid2, float* wt2,
 }
 
 void becke_weight_2d(int gs, double* grid1, double* wt1, double* grid2, double* wt2,
-                     double zeta1, double zeta2, double A2, double B2, double C2)
+                     bool is_Z, double zeta1, double zeta2, double A2, double B2, double C2)
 {
   int mu_order = 3;
 
@@ -962,8 +962,19 @@ void becke_weight_2d(int gs, double* grid1, double* wt1, double* grid2, double*
   double R  = sqrt(A2*A2+B2*B2+C2*C2);
   const double oR = 1./R;
 
-  const double a1 = becke_a_zeta(zeta1,zeta2);
-  const double a2 = becke_a_zeta(zeta2,zeta1);
+  double a1 = 0.;
+  double a2 = 0.;
+  if (is_Z)
+  {
+   //not actually zeta input, but Z input to this ftn
+    a1 = becke_a(zeta1,zeta2);
+    a2 = becke_a(zeta2,zeta1);
+  }
+  else
+  {
+    a1 = becke_a_zeta(zeta1,zeta2);
+    a2 = becke_a_zeta(zeta2,zeta1);
+  }
 
   //printf(" a1/2: %8.5f %8.5f \n",a1,a2);
 
@@ -5220,9 +5231,6 @@ void compute_delt(int natoms, int* atno, double* coords, bool gbasis, vector<vec
 
   if (sgs_basis) { printf("\n ERROR: compute_rho for SGS basis does not support drho \n"); exit(-1); }
 
-  if (gbasis)
-    printf("  TESTING: compute_delt for gbasis \n");
-
   int N = basis.size();
   int N2 = N*N;
   int* n2i = new int[natoms];
diff --git a/src/integrals/cpu_util.cpp b/src/integrals/cpu_util.cpp
index e14600d..3f44e36 100644
--- a/src/integrals/cpu_util.cpp
+++ b/src/integrals/cpu_util.cpp
@@ -34,6 +34,19 @@ extern void dsyevx_(
     int* iwork, int* IFAIL,
     int* info );
 
+extern void dsygvx_(
+    int* itype, char* jobz, char* range, char* uplo,
+    int* n,
+    double* A, int* lda,
+    double* B, int* ldb,
+    double* vl, double* vu, int* il, int* iu,
+    double* abstol, int* m,
+    double* W,
+    double* Z, int* ldz,
+    double* work, int* lwork,
+    int* iwork, int* IFAIL,
+    int* info );
+
 extern void dgetrf_(int* M, int* N, double* A, int* LDA, int* IPIV, int* INFO);
 
 extern void dgetri_(int* N, double* A, int* LDA, int* IPIV, double* WORK, int* LWORK, int* INFO);
@@ -44,6 +57,11 @@ extern "C" {
                 double* b, int* ldb, int* info);
 }
 
+extern "C" {
+    void dgels_(char* trans, int* m, int* n, int* nrhs, double* a, int* lda,
+                double* b, int* ldb, double* work, int* lwork, int* info);
+}
+
 void print_square(int N, double* A);
 void print_square_sm(int N, double* A);
 
@@ -58,6 +76,31 @@ double randomf(double a, double b)
   return randn;
 }
 
+void solve_axeb_ls(int m, int n, double* A, double* b)
+{
+  //int m = rows;    // number of rows
+  //int n = cols;    // number of columns
+  int nrhs = 1;
+  int lwork = -1;
+  int info;
+
+ // A is m*n (column-major), b is max(m,n)
+ // b must be allocated to max(m,n) to hold the solution in-place
+  int size = (int){m>n ? m:n};
+  //double b[m > n ? m : n];
+
+  double work_query;
+  dgels_("T", &m, &n, &nrhs, A, &m, b, &size, &work_query, &lwork, &info);
+
+  lwork = (int)work_query;
+  double* work = new double[lwork];
+  dgels_("T", &m, &n, &nrhs, A, &m, b, &size, work, &lwork, &info);
+
+  delete [] work;
+
+  return;
+}
+
 void solve_axeb(int dim, double* A, double* b)
 {
   int info;
@@ -68,6 +111,35 @@ void solve_axeb(int dim, double* A, double* b)
   return;
 }
 
+double determinant(int N, double* A)
+{
+  int* ipiv = new int[N];
+  int LDA = N;
+  int info = 0;
+
+  dgetrf_(&N,&N,A,&LDA,ipiv,&info);
+
+  if (info<0)
+    printf(" dgetrf failed \n");
+
+  double det = 1.0;
+  if (info==0)
+  {
+    for (int i=0;i<N;i++)
+    {
+      det *= A[i*N+i];
+      if (i+1<N && ipiv[i]!=i+1)
+        det *= -1.0; // row swap flips sign
+    }
+  }
+  else
+    det = 0.; //singular
+
+  delete [] ipiv;
+
+  return det;
+}
+
 void expmat_complex_cpu(int N, double* theta, double* thetai, double* etheta)
 {
   printf("\n ERROR expmat_complex_cpu broken by new compiler \n"); exit(-1);
@@ -270,6 +342,63 @@ int la_diagR(int size, double* A, double* eigen, double* eigeni)
   return 0;
 }
 
+//overwrites A+S
+void la_diag_gen(int neig, int s1, double* A, double* S, double* Ae)
+{
+  int s1a = s1;
+  int LDA = s1;
+  int LDB = s1;
+
+  int itype = 1;  // A*x = lambda*S*x
+
+  char jobz = 'V';
+  char range = 'I';
+  if (neig==s1)
+    range = 'A';
+  char uplo = 'u';
+
+  int il = 1;
+  int iu = neig;
+
+ //abstol not necessarily zero. see LAPACK manual
+  double abstol = 0.;
+  double vl = 1.; double vu = -1.;
+  int neval = s1;
+
+  double* EVec = new double[LDA*s1];
+
+  int LenWork = 32*s1; //LWORK >= max(1,8*N)
+  double* Work = new double[LenWork];
+
+  int LenIWork = 10*s1; //IWORK needs 5*N min
+  int* IWork = new int[LenIWork];
+  int* IFail = new int[s1];
+
+  int Info = 0;
+
+  dsygvx_(&itype,&jobz,&range,&uplo,&s1a,A,&LDA,S,&LDB,&vl,&vu,&il,&iu,&abstol,
+       &neval,Ae,EVec,&LDA,Work,&LenWork,IWork,IFail,&Info);
+
+  if (Info!=0)
+  {
+    if (Info > s1)
+      printf(" ERROR: S matrix not positive definite (leading minor %i), Info = %i \n",Info-s1,Info);
+    else
+      printf(" la_diag_gen Info = %d\n",Info);
+  }
+
+  for (int i=0;i<s1;i++)
+  for (int j=0;j<s1;j++)
+    A[i*s1+j] = EVec[i*s1+j];
+
+  delete [] EVec;
+  delete [] Work;
+  delete [] IWork;
+  delete [] IFail;
+
+  return;
+}
+
 void la_diag(int neig, int s1, double* A, double* Ae)
 {
   #define USE_DSYEVX 1
@@ -347,6 +476,7 @@ void la_diag(int neig, int s1, double* A, double* Ae)
   delete [] IWork;
   delete [] IFail;
 
+  return;
 }
 
 void test_diag()
@@ -1101,10 +1231,10 @@ void mat_times_mat_bt_cpu(double* C, double* A, double* B, int N)
   int LDA = N;
   int LDB = N;
   int LDC = N;
-   
+
   double ALPHA = 1.0;
   double BETA = 0.0;
-  
+
  //C := alpha*op( A )*op( B ) + beta*C (op means A or B, possibly transposed)
  //CBlas version
   //cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasTrans,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC);
diff --git a/src/integrals/cuda_util.cpp b/src/integrals/cuda_util.cpp
index 1e8b47f..fbc906d 100644
--- a/src/integrals/cuda_util.cpp
+++ b/src/integrals/cuda_util.cpp
@@ -494,7 +494,7 @@ void trans(double* Bt, double* B, int m, int n)
   return;
 }
 
-int mat_root_inv_stable_cusolver(double* A, int size, double delta, cusolverDnHandle_t& cu_hdl)
+int mat_root_inv_stable_cusolver(double* A, int size, double delta, cusolverDnHandle_t& cu_hdl, int prl)
 {
   int s2 = size*size;
 
@@ -513,28 +513,32 @@ int mat_root_inv_stable_cusolver(double* A, int size, double delta, cusolverDnHa
 
   diagonalize_cusolver(size,size,B,Beigen,cu_hdl);
 
- #if 1
-  #pragma acc update self(Beigen[0:size])
-  printf(" mat_root_inv_stable eigenvalues:");
-  for (int i=0;i<std::min(10,size);i++)
-  if (fabs(Beigen[i])<1e10)
-    printf(" %5.2e",Beigen[i]);
-  else
-    printf(" large");
-  printf("\n");
- #endif
-
   int nsmall = 0;
  #pragma acc parallel loop present(Beigen[0:size]) reduction(+:nsmall)
   for (int i=0;i<size;i++)
   if (Beigen[i]<delta)
-  {
-    Beigen[i] = 1.e20;
     nsmall++;
+
+  if (prl>1)
+  {
+    #pragma acc update self(Beigen[0:size])
+    printf(" mat_root_inv_stable eigenvalues:");
+    for (int i=0;i<std::max(nsmall,std::min(10,size));i++)
+    if (fabs(Beigen[i])<1e10)
+      printf(" %5.2e",Beigen[i]);
+    else
+      printf(" large");
+    printf("\n");
   }
-  printf("  found %i small eigen \n",nsmall);
 
-  if (nsmall>0)
+ #pragma acc parallel loop present(Beigen[0:size])
+  for (int i=0;i<nsmall;i++)
+    Beigen[i] = 1.e40;
+
+  if (prl>0)
+    printf("  found %i small eigen \n",nsmall);
+
+  if (prl>0 && nsmall>0)
   {
     printf(" lowest vector:");
     #pragma acc update self(B[0:s2])
@@ -580,6 +584,11 @@ int mat_root_inv_stable_cusolver(double* A, int size, double delta, cusolverDnHa
   return nsmall;
 }
 
+int mat_root_inv_stable_cusolver(double* A, int size, double delta, cusolverDnHandle_t& cu_hdl)
+{
+  return mat_root_inv_stable_cusolver(A,size,delta,cu_hdl,1);
+}
+
 int mat_root_inv_cusolver(double* A, int size, cusolverDnHandle_t& cu_hdl)
 {
   int s2 = size*size;
diff --git a/src/integrals/gauss.cpp b/src/integrals/gauss.cpp
index 040f317..eb61072 100644
--- a/src/integrals/gauss.cpp
+++ b/src/integrals/gauss.cpp
@@ -9,490 +9,181 @@
 //eval_gh: check factor of sqrt(2)?
 
 
-void eval_vghd(int gsa, double* grid, double* val, int n, int l, int m, double norm, double zt)
+void eval_vghd(int gsa, double* grid, double* val, int n, int l, double norm, double zt)
 {
   int gsa6 = 6*gsa;
-  double nph = l+0.5;
 
   if (l>5) { printf("\n ERROR: l>6 in eval_vgh \n"); exit(-1); }
 
- #if 1
-  double normf = 2.*PI/(2.*l+1.)*norm;
-  double norm1 = pow(zt,-0.5*(2.*l+3.));
-  double norm2 = pow(zt,-1.-0.5*l);
-  double gf1 = 0.5*(2.*l+3.);
-  double gf2 = 1.-0.5*l;
-
-  #pragma acc parallel loop present(grid[0:gsa6],val[0:gsa])
-  for (int j=0;j<gsa;j++)
-  {
-    double r = grid[6*j+3];
-    double zr2 = zt*r*r;
-    double rl1 = pow(r,-l-1.);
-    double rl  = pow(r,l);
-    val[j] = normf*(norm1*rl1*igam(gf1,zr2) + norm2*rl*igamc(gf2,zr2));
-  }
- #endif
-
- #if 0
-  double normf = 0.5*PI*norm/zt;
-  //normf *= 11.136655993664;
- //CPMZ check all of these
- /*
-  if (l==1) normf *= 0.5;
-  if (l==2) normf *= 0.75;
-  if (l==3) normf *= 1.875;
-  if (l==4) normf *= 6.5625;
-  if (l==5) normf *= 29.53125;
-  if (l==6) normf *= 162.421875;
- */
+  double nph = l+0.5;
+  double normf = 0.25*PI*PI*norm/zt;
+  //double normf = 0.5*norm/zt; //this norm gives A matrix same as eri_2c
+  normf *= tgamma(nph);
 
  #pragma acc parallel loop present(grid[0:gsa6],val[0:gsa])
   for (int j=0;j<gsa;j++)
   {
     double r = grid[6*j+3];
     double zr2 = zt*r*r;
-    //double rl = pow(r,l);
 
     double Fn = pow(zr2,-nph)*igam(nph,zr2);
     val[j] = normf*Fn;
   }
- #endif
 
   return;
 }
 
-#if 0
-void integrate_hole_para_gh(double* rdm, bool full_rdm, bool hfx_on, int Nc, int No, int M, int natoms, int* atno, double* coords, int gs, int gsb, vector<vector<double> >& basis,
-                    double* Pao, double* Pmo, double* jCA, float* grid, float* gridb, float* wt, float* wtb, float* rho, float* vxch, int prl)
+//CPMZ update this to work with MO basis
+void gen_4cj_overlap(double* ol4j, int nrad, int nang, double* ang_g, double* ang_w,
+                     int natoms, int* atno, double* coords, vector<vector<double> > basis, int prl)
 {
-  #define OUTER_DIVIDE 1
-
-  int ngpu = acc_get_num_devices(acc_device_nvidia);
-  if (ngpu<1) ngpu = 1;
+  if (prl>-1) printf(" TESTING: gen_4cj_overlap \n");
 
-  printf("\n creating xc hole potential in Gaussian basis (ngpu: %i) \n",ngpu);
-  fflush(stdout);
-
-  if (Nc>0)
-  {
-    printf("\n ERROR: Nc must be zero \n");
-    exit(1);
-  }
+  int gs = nrad*nang;
+  int gsa = gs*natoms;
+  int gsa3 = gsa*3;
+  int gsa6 = gsa*6;
 
   int N = basis.size();
   int N2 = N*N;
-
-  int M2 = M*M;
-  int M3 = M*M2;
-  int M4 = M2*M2;
-
-  float* Pmof = new float[N2];
-  for (int m=0;m<N2;m++)
-    Pmof[m] = Pmo[m];
-
-  if (M<8 && prl>-2)
-  {
-    printf("\n jCA: \n");
-    print_square(N,jCA);
-    printf("\n Pao: \n");
-    print_square(N,Pao);
-    printf("\n Pmo: \n");
-    print_square(N,Pmo);
-    print_rdm(M,rdm);
-  }
-
-  if (!full_rdm)
-  {
-   //removes one operation from within loop
-    for (int p=0;p<M;p++)
-    for (int q=0;q<M;q++)
-    for (int r=0;r<M;r++)
-    for (int s=0;s<M;s++)
-      rdm[p*M3+q*M2+r*M+s] -= Pmo[p*N+q]*Pmo[r*N+s];
-  }
-
-  float* rdmf = new float[M4];
-  for (int m=0;m<M4;m++)
-    rdmf[m] = rdm[m];
-
-  int gsa = natoms*gs;
-  int gsa6 = 6*gsa;
-  int gsba = natoms*gsb;
-  int gsba6 = 6*gsba;
+  int N3 = N2*N;
+  int N4 = N3*N;
 
   int* n2i = new int[natoms];
   int imaxN = get_imax_n2i(natoms,N,basis,n2i);
 
   int iN = N;
-
- //tmp arrays for summing over contracted Gaussians
-  float* val0 = new float[gsa];
-  float* val0b = new float[gsba];
-
- // 1 --> AO
- // 2 --> MO
-  float** val1 = new float*[iN];
-  for (int i=0;i<iN;i++)
-    val1[i] = new float[gsa];
-  float** val2 = new float*[iN];
-  for (int i=0;i<iN;i++)
-    val2[i] = new float[gsa];
-
-  float** val1b = new float*[iN];
+  double* val0 = new double[gsa];
+  double* val0p = new double[gsa3];
+  double** val1  = new double*[iN];
+  double** val1p = new double*[iN];
   for (int i=0;i<iN;i++)
-    val1b[i] = new float[gsba];
-  float** val2b = new float*[iN];
+    val1[i]  = new double[gsa];
   for (int i=0;i<iN;i++)
-    val2b[i] = new float[gsba];
-
-  float* grid1 = new float[gsa6];
-  float* grid1b = new float[gsba6];
-
-  float* vxcht = new float[gsa];
-  float* rdm2 = new float[M2];
-
-  #pragma omp parallel for schedule(static) num_threads(ngpu)
-  for (int g=0;g<ngpu;g++)
-  {
-    acc_set_device_num(g,acc_device_nvidia);
+    val1p[i] = new double[gsa3];
 
-    if (g>0)
-    {
-      #pragma acc enter data create(grid[0:gsa6])
-    }
-    #pragma acc enter data create(vxcht[0:gsa])
+  double* grid = new double[gsa6];
+  double* wt = new double[gsa];
+  double* grid1 = new double[gsa6];
+  #pragma acc enter data create(grid[0:gsa6],wt[0:gsa],grid1[0:gsa6])
 
-    #pragma acc enter data copyin(rdmf[0:M4],Pmof[0:N2])
-    #pragma acc enter data create(rdm2[0:M2])
+  #pragma acc enter data copyin(ang_g[0:3*nang],ang_w[0:nang])
+  get_becke_grid_full(natoms,atno,coords,nrad,nang,ang_g,ang_w,6,grid,wt);
+  #pragma acc exit data delete(ang_g[0:3*nang],ang_w[0:nang])
 
-    #pragma acc enter data copyin(jCA[0:N2],Pmo[0:N2])
-    #pragma acc enter data create(val0[0:gsa],val0b[0:gsba],val1[0:iN][0:gsa],val2[0:iN][0:gsa],val1b[0:iN][0:gsba],val2b[0:iN][0:gsba])
+  #pragma acc enter data create(val0[0:gsa],val1[0:iN][0:gsa])
+  #pragma acc enter data create(val0p[0:gsa3],val1p[0:iN][0:gsa3])
 
-    if (g==0)
-    {
-      #pragma acc enter data create(grid1[0:gsa6],grid1b[0:gsba6])
-    }
+  for (size_t j=0;j<N2*N2;j++)
+    ol4j[j] = 0.;
 
-    #pragma acc parallel loop present(vxcht[0:gsa])
-    for (int j=0;j<gsa;j++)
-      vxcht[j] = 0.;
-  }
-  acc_set_device_num(0,acc_device_nvidia);
+  #pragma acc parallel loop collapse(2) present(val1[0:iN][0:gsa])
+  for (int i1=0;i1<N;i1++)
+  for (int j=0;j<gsa;j++)
+    val1[i1][j] = 0.;
 
+  #pragma acc parallel loop collapse(2) present(val1p[0:iN][0:gsa3])
+  for (int i1=0;i1<N;i1++)
+  for (int j=0;j<gsa3;j++)
+    val1p[i1][j] = 0.;
 
-  const int ig = 10; //first index to exp in Gaussian basis
+  const int ig = 10;
 
   for (int m=0;m<natoms;m++)
   {
    //working on this block of the matrix
     int s1 = 0; if (m>0) s1 = n2i[m-1]; int s2 = n2i[m];
 
-    float Z1 = (float)atno[m];
     double A1 = coords[3*m+0]; double B1 = coords[3*m+1]; double C1 = coords[3*m+2];
 
     copy_grid(gsa,grid1,grid);
     recenter_grid_zero(gsa,grid1,-A1,-B1,-C1);
-    copy_grid(gsba,grid1b,gridb);
-    recenter_grid_zero(gsba,grid1b,-A1,-B1,-C1);
-
-    #pragma acc parallel loop collapse(2) present(val1[0:iN][0:gsa])
-    for (int i1=s1;i1<s2;i1++)
-    for (int j=0;j<gsa;j++)
-      val1[i1][j] = 0.f;
-
-    #pragma acc parallel loop collapse(2) present(val1b[0:iN][0:gsba])
-    for (int i1=s1;i1<s2;i1++)
-    for (int j=0;j<gsba;j++)
-      val1b[i1][j] = 0.f;
 
     for (int i1=s1;i1<s2;i1++)
     {
-      int ii1 = i1;
-
       vector<double> basis1 = basis[i1];
       int l1 = basis1[1]; int m1 = basis1[2]; int ng = basis1[3];
       int in = ig + ng; //index to find norm
 
-      float* valm = val1[ii1];
-      float* valn = val1b[ii1];
+      double* valm = val1[i1];
+      double* valmp = val1p[i1];
       for (int j=0;j<ng;j++)
       {
         double zeta1 = basis1[ig+j]; double norm1 = basis1[in+j];
 
-        eval_gh(gsa,grid1,val0,l1,m1,norm1,zeta1);
+        eval_ghd(gsa,grid1,val0,l1,m1,norm1,zeta1);
 
-        #pragma acc parallel loop present(val0[0:gsa],valm[0:gsa])
+       #pragma acc parallel loop present(val0[0:gsa],valm[0:gsa])
         for (int k=0;k<gsa;k++)
           valm[k] += val0[k];
 
-        eval_gh(gsba,grid1b,val0b,l1,m1,norm1,zeta1);
+        eval_pd_gh(gsa,grid1,val0p,l1+1,l1,m1,norm1,zeta1);
 
-        #pragma acc parallel loop present(val0b[0:gsba],valn[0:gsba])
-        for (int k=0;k<gsba;k++)
-          valn[k] += val0b[k];
+       #pragma acc parallel loop present(val0p[0:gsa3],valmp[0:gsa3])
+        for (int k=0;k<gsa3;k++)
+          valmp[k] += val0p[k];
       }
-    }
-  }
-
 
- #if OUTER_DIVIDE
- //need rho
-  #pragma acc parallel loop present(rho[0:gsa])
-  for (int j=0;j<gsa;j++)
-    rho[j] = RHO_ZERO;
+    } //loop i1 over basis
+  } //loop m over natoms
 
+ //compute all
   for (int i1=0;i1<N;i1++)
+  for (int i2=0;i2<N;i2++)
   {
-    int ii1 = i1;
-    float* valn = val1[ii1];
+    double* valm = val1[i1];
+    double* valnp = val1p[i2];
+    if (i2==0 && N>20 && i1%2) printf(" . ");
 
-    for (int i2=0;i2<N;i2++)
+    for (int i3=0;i3<N;i3++)
+    for (int i4=0;i4<N;i4++)
     {
-      float d1 = Pao[i1*N+i2];
+      int ind = i1*N3+i2*N2+i3*N+i4;
 
-      int ii2 = i2;
-      float* valm = val1[ii2];
+      double* valp = val1[i3];
+      double* valqp = val1p[i4];
 
-      #pragma acc parallel loop present(rho[0:gsa],valm[0:gsa],valn[0:gsa])
+      double vtx = 0.; double vty = 0.; double vtz = 0.;
+     #pragma acc parallel loop present(valm[0:gsa],valnp[0:gsa3],valp[0:gsa],valqp[0:gsa3]) reduction(+:vtx,vty,vtz)
       for (int j=0;j<gsa;j++)
-        rho[j] += d1*valn[j]*valm[j];
-    }
-  }
- #endif
-
- //AO-->MO
-  #pragma acc parallel loop collapse(2) present(val2[0:iN][0:gsa])
-  for (int i1=0;i1<N;i1++)
-  for (int j=0;j<gsa;j++)
-    val2[i1][j] = 0.f;
-  #pragma acc parallel loop collapse(2) present(val2b[0:iN][0:gsba])
-  for (int i1=0;i1<N;i1++)
-  for (int j=0;j<gsba;j++)
-    val2b[i1][j] = 0.f;
-
-  #pragma acc parallel loop collapse(2) present(val1[0:iN][0:gsa],val2[0:iN][0:gsa],jCA[0:N2])
-  for (int i1=0;i1<N;i1++)
-  for (int j=0;j<gsa;j++)
-  {
-    float v2 = 0.f;
-    #pragma acc loop reduction(+:v2)
-    for (int k=0;k<N;k++)
-      v2 += jCA[k*N+i1]*val1[k][j];
-    val2[i1][j] = v2;
-  }
-  #pragma acc parallel loop collapse(2) present(val1b[0:iN][0:gsba],val2b[0:iN][0:gsba],jCA[0:N2])
-  for (int i1=0;i1<N;i1++)
-  for (int j=0;j<gsba;j++)
-  {
-    float v2 = 0.f;
-    #pragma acc loop reduction(+:v2)
-    for (int k=0;k<N;k++)
-      v2 += jCA[k*N+i1]*val1b[k][j];
-    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;
-  double den = 0.;
-  #pragma omp parallel for schedule(static) num_threads(ngpu) reduction(+:den)
-  for (int g=0;g<ngpu;g++)
-  {
-    acc_set_device_num(g,acc_device_nvidia);
-    int s1 = g*gdiv;
-    int s2 = (g+1)*gdiv; if (g==ngpu-1) s2 = gsba;
-
-    if (prl>1) printf("  s1/2: %6i %6i \n",s1,s2);
-    for (int j=s1;j<s2;j++)
-    {
-      float x1 = gridb[6*j+0]; float y1 = gridb[6*j+1]; float z1 = gridb[6*j+2];
-      float wt1 = wtb[j];
-
-      #pragma acc parallel loop present(rdm2[0:M2])
-      for (int m=0;m<M2;m++)
-        rdm2[m] = 0.f;
-
-      float dtt = 0.f;
-      #pragma acc parallel loop collapse(2) present(rdm2[0:M2],rdmf[0:M4],Pmof[0:N2],val2b[0:iN][0:gsba]) reduction(+:dtt)
-      for (int r=0;r<M;r++)
-      for (int s=0;s<M;s++)
-      {
-        float dt = 0.f;
-        //float p1 = Pmof[r*N+s];
-
-        #pragma acc loop collapse(2) reduction(+:dt)
-        for (int p=0;p<M;p++)
-        for (int q=0;q<M;q++)
-        {
-          float d0 = rdmf[p*M3+q*M2+r*M+s];
-          d0 *= val2b[p][j]*val2b[q][j];
-          dt += d0;
-        }
-        rdm2[r*M+s] = dt;
-        dtt += dt;
-      }
-
-      double d1b = 0.;
-      #pragma acc parallel loop collapse(2) present(Pmof[0:N2],val2b[0:iN][0:gsba]) reduction(+:d1b)
-      for (int p=0;p<N;p++)
-      for (int q=0;q<N;q++)
       {
-        float d0 = Pmof[p*N+q];
-        d0 *= val2b[p][j]*val2b[q][j];
+        double v10 = valm[j]*valp[j];
+        double v1x = valnp[3*j+0]*valqp[3*j+0];
+        double v1y = valnp[3*j+1]*valqp[3*j+1];
+        double v1z = valnp[3*j+2]*valqp[3*j+2];
 
-        d1b += d0;
+        vtx += v10*v1x;
+        vty += v10*v1y;
+        vtz += v10*v1z;
       }
-      den += d1b*wt1;
-
-      //if (fabs(dtt)>1.e-8 && fabs(d1b)>1.e-8)
-      #pragma acc parallel loop present(vxcht[0:gsa],Pmo[0:N2],Pmof[0:N2],rdm2[0:M2],grid[0:gsa6],val2[0:iN][0:gsa],val2b[0:iN][0:gsba])
-      for (int k=0;k<gsa;k++)
-      {
-        float x12 = grid[6*k+0]-x1; float y12 = grid[6*k+1]-y1; float z12 = grid[6*k+2]-z1;
-
-        float rn = sqrtf(x12*x12+y12*y12+z12*z12)+RDEN;
-        float or1 = 1.f/rn;
-
-       #if !OUTER_DIVIDE
-        double d1 = RHO_ZERO;
-        #pragma acc loop collapse(2) reduction(+:d1)
-        for (int p=0;p<N;p++)
-        for (int q=0;q<N;q++)
-        {
-          float d0 = Pmof[p*N+q];
-          d0 *= val2[p][k]*val2[q][k];
-
-          d1 += d0;
-        }
-       #endif
-
-        double d2 = 0.;
-        //if (fabs(d1)>1.e-8)
-        #pragma acc loop collapse(2) reduction(+:d2)
-        for (int r=0;r<M;r++)
-        for (int s=0;s<M;s++)
-        {
-          float d0 = rdm2[r*M+s];
-          d0 *= val2[r][k]*val2[s][k];
-          d2 += d0;
-        }
-       #if !OUTER_DIVIDE
-        d2 /= d1;
-       #endif
 
-        vxcht[k] += d2*wt1*or1;
-
-      } //loop k (on gpu)
-    } //loop j
-
-  } //loop g (omp)
-
-  for (int j=0;j<gsa;j++)
-    vxch[j] = 0.;
-  //int gdiv2 = gsa/ngpu+1;
-  for (int g=0;g<ngpu;g++)
-  {
-    acc_set_device_num(g,acc_device_nvidia);
-    #pragma acc update self(vxcht[0:gsa])
-
-    //#pragma omp parallel for schedule(static,gdiv2) num_threads(ngpu)
-    for (int k=0;k<gsa;k++)
-      vxch[k] += vxcht[k];
-  }
-  acc_set_device_num(0,acc_device_nvidia);
-  #pragma acc update device(vxch[0:gsa])
-
-  double Pao_den = 0.;
- #if OUTER_DIVIDE
-  #pragma acc parallel loop present(rho[0:gsa],wt[0:gsa]) reduction(+:Pao_den)
-  for (int k=0;k<gsa;k++)
-    Pao_den += rho[k]*wt[k];
-
-  #pragma acc parallel loop present(vxch[0:gsa],rho[0:gsa])
-  for (int k=0;k<gsa;k++)
-    vxch[k] /= rho[k];
-
-  #pragma acc update self(vxch[0:gsa])
- #endif
-
-
-  printf("\n");
-  printf("  total den (Pmo/Pao): %9.6f %9.6f \n",den,Pao_den);
-
-
-  if (full_rdm)
-    printf("\n vxch(full): \n");
-  else
-    printf("\n vxch: \n");
-  for (int n=0;n<natoms;n++)
-  {
-    for (int m=0;m<gs;m+=pr_inc)
-      printf(" %8.5f",vxch[n*gs+m]);
-    printf("\n");
-  }
-  printf("\n");
-
- //cleanup
-  #pragma omp parallel for schedule(static) num_threads(ngpu)
-  for (int g=0;g<ngpu;g++)
-  {
-    acc_set_device_num(g,acc_device_nvidia);
-
-    if (g>0)
-    {
-      #pragma acc exit data delete(grid[0:gsa6])
-    }
-    #pragma acc exit data delete(vxcht[0:gsa])
-    #pragma acc exit data delete(rdmf[0:M4])
-    #pragma acc exit data delete(rdm2[0:M2],Pmof[0:N2])
-    #pragma acc exit data delete(jCA[0:N2],Pmo[0:N2])
-    #pragma acc exit data delete(val0[0:gsa],val0b[0:gsba],val1[0:iN][0:gsa],val2[0:iN][0:gsa],val1b[0:iN][0:gsba],val2b[0:iN][0:gsba])
-    if (g==0)
-    {
-      #pragma acc exit data delete(grid1[0:gsa6],grid1b[0:gsba6])
+      ol4j[ind] = vtx;
+      ol4j[N4+ind] = vty;
+      ol4j[2*N4+ind] = vtz;
     }
   }
 
-  delete [] n2i;
+  #pragma acc exit data delete(grid[0:gsa6],wt[0:gsa],grid1[0:gsa6])
+  #pragma acc exit data delete(val0[0:gsa],val1[0:iN][0:gsa])
+  #pragma acc exit data delete(val0p[0:gsa3],val1p[0:iN][0:gsa3])
 
+  delete [] grid;
+  delete [] wt;
   delete [] grid1;
-  delete [] grid1b;
-
-  delete [] val0;
-  delete [] val0b;
+  delete [] n2i;
 
   for (int i=0;i<iN;i++)
     delete [] val1[i];
   delete [] val1;
-
-  for (int i=0;i<iN;i++)
-    delete [] val2[i];
-  delete [] val2;
-
   for (int i=0;i<iN;i++)
-    delete [] val1b[i];
-  delete [] val1b;
-
-  for (int i=0;i<iN;i++)
-    delete [] val2b[i];
-  delete [] val2b;
-
-  delete [] Pmof;
-  delete [] vxcht;
+    delete [] val1p[i];
+  delete [] val1p;
 
-  delete [] rdmf;
+  delete [] val0;
+  delete [] val0p;
 
   return;
 }
-#endif
 
 void wf_to_grid_gh(bool divide_by_rho, bool calc_rho, int natoms, int* atno, double* coords, vector<vector<double> > basis,
                    double* Pao, double* gfao, int gs, double* grid, double* wt,
@@ -561,10 +252,10 @@ void wf_to_grid_gh(bool divide_by_rho, bool calc_rho, int natoms, int* atno, dou
   double* grid2 = new double[gsa6];
 
   #pragma acc enter data create(grid1[0:gsa6],grid2[0:gsa6])
-  #pragma acc enter data create(val0[0:gsa],val1[0:iN][0:gsa],val2[0:iN][0:gsa])
+  #pragma acc enter data create(val0[0:gsa],val0p[0:gsa3],val1[0:iN][0:gsa],val2[0:iN][0:gsa])
   if (need_g)
   {
-    #pragma acc enter data create(val0p[0:gsa3],val1p[0:iN][0:gsa3],val2p[0:iN][0:gsa3])
+    #pragma acc enter data create(val1p[0:iN][0:gsa3],val2p[0:iN][0:gsa3])
   }
   if (need_L)
   {
@@ -820,6 +511,7 @@ void wf_to_grid_gh(bool divide_by_rho, bool calc_rho, int natoms, int* atno, dou
   }
 
   if (prl>0 || calc_rho)
+  if (rho!=NULL)
   {
     double dtot = 0.;
    #pragma acc parallel loop present(rho[0:gsa],wt[0:gsa])
@@ -879,10 +571,10 @@ void wf_to_grid_gh(bool divide_by_rho, bool calc_rho, int natoms, int* atno, dou
     gfao[j] *= 0.5;
 
   #pragma acc exit data delete(grid1[0:gsa6],grid2[0:gsa6])
-  #pragma acc exit data delete(val0[0:gsa],val1[0:iN][0:gsa],val2[0:iN][0:gsa])
+  #pragma acc exit data delete(val0[0:gsa],val0p[0:gsa3],val1[0:iN][0:gsa],val2[0:iN][0:gsa])
   if (need_g)
   {
-    #pragma acc exit data delete(val0p[0:gsa3],val1p[0:iN][0:gsa3],val2p[0:iN][0:gsa3])
+    #pragma acc exit data delete(val1p[0:iN][0:gsa3],val2p[0:iN][0:gsa3])
   }
   if (need_L)
   {
diff --git a/src/integrals/integrals.cpp b/src/integrals/integrals.cpp
index ffaf609..285a514 100644
--- a/src/integrals/integrals.cpp
+++ b/src/integrals/integrals.cpp
@@ -4242,7 +4242,7 @@ void reduce_Exyz(int i1, int i2, int N, int gs, double* val1m, double* val2m, do
   int gs6 = 6*gs;
   int N2 = N*N;
 
-  const double a = 11.; //33, 66 also viable?
+  const double a = 25.; //33, 66 also viable?
 
   double valx = 0.; double valy = 0.; double valz = 0.;
   #pragma acc parallel loop present(val1m[0:gs],val2m[0:gs],grid1m[0:gs6]) reduction(+:valx,valy,valz)
@@ -4281,7 +4281,7 @@ void reduce_Exyz_2(int i1, int i2, int N, int gs, double* val1m, double* val1n,
   int gs6 = 6*gs;
   int N2 = N*N;
 
-  const double a = 11.;
+  const double a = 25.;
 
   double valx = 0.; double valy = 0.; double valz = 0.;
   #pragma acc parallel loop present(val1m[0:gs],val1n[0:gs],val2m[0:gs],val2n[0:gs],grid1m[0:gs6],grid2n[0:gs6]) reduction(+:valx,valy,valz)
@@ -4357,7 +4357,7 @@ void compute_Exyz(double rconf, double pconf, int natoms, int* atno, double* coo
 
   int* n2i = new int[natoms];
   int imaxN = get_imax_n2i(natoms,N,basis,n2i);
-  printf("  iN: %i \n",imaxN);
+  //printf("  iN: %i \n",imaxN);
 
   const int ig = 10;
 
@@ -4484,6 +4484,7 @@ void compute_Exyz(double rconf, double pconf, int natoms, int* atno, double* coo
    //if (natoms>1) { printf("  WARNING: testing >1 atom in compute_Exyz \n"); }
 
    //complete but needs testing
+    if (1)
     for (int n=m+1;n<natoms;n++)
     {
       int s3 = 0; if (n>0) s3 = n2i[n-1]; int s4 = n2i[n];
@@ -4505,7 +4506,7 @@ void compute_Exyz(double rconf, double pconf, int natoms, int* atno, double* coo
 
        //optimize this
         //becke_weight_2d(gs,grid1m,wt1,grid2m,wt2,zeta1,zeta2,A12,B12,C12);
-        becke_weight_2d(gs,grid1m,wt1,grid2m,wt2,Z1,Z2,A12,B12,C12);
+        becke_weight_2d(gs,grid1m,wt1,grid2m,wt2,1,Z1,Z2,A12,B12,C12);
 
         copy_grid(gs,grid2n,grid2m);
         recenter_grid(gs,grid2n,-A12,-B12,-C12);      //grid 2 centered on atom 1
@@ -4517,7 +4518,7 @@ void compute_Exyz(double rconf, double pconf, int natoms, int* atno, double* coo
         add_r1_to_grid(gs,grid2m,0.,0.,0.);
 
         for (int i1=s1;i1<s2;i1++)
-        for (int i2=s1;i2<=i1;i2++)
+        for (int i2=s3;i2<s4;i2++)
         {
           vector<double> basis1 = basis[i1];
           int n1 = basis1[0]; int l1 = basis1[1]; int m1 = basis1[2]; int ng1 = basis1[3];
@@ -4595,7 +4596,7 @@ void compute_Exyz(double rconf, double pconf, int natoms, int* atno, double* coo
         recenter_grid(gs,grid2m,A12,B12,C12);
 
        //optimize this
-        becke_weight_2d(gs,grid1m,wt1,grid2m,wt2,zeta1,zeta2,A12,B12,C12);
+        becke_weight_2d(gs,grid1m,wt1,grid2m,wt2,0,zeta1,zeta2,A12,B12,C12);
         //becke_weight_2d(gs,grid1m,wt1,grid2m,wt2,Z1,Z2,A12,B12,C12);
 
         copy_grid(gs,grid2n,grid2m);
@@ -4865,7 +4866,7 @@ void compute_Sd(int natoms, int* atno, float* coords, vector<vector<double> > &b
         recenter_grid(gs,grid2m,A12,B12,C12);
 
        //optimize this
-        becke_weight_2d(gs,grid1m,wt1,grid2m,wt2,zeta1,zeta2,A12,B12,C12);
+        becke_weight_2d(gs,grid1m,wt1,grid2m,wt2,0,zeta1,zeta2,A12,B12,C12);
         //becke_weight_2d(gs,grid1m,wt1,grid2m,wt2,Z1,Z2,A12,B12,C12);
 
         copy_grid(gs,grid2n,grid2m);
@@ -5544,7 +5545,7 @@ void compute_all_2c_v2d(bool do_overlap, int natoms, int* atno, float* coords, v
 
       acc_copy(gs,wtt1,wt1);
       //becke_weight_2c(gs,grid1,wtt1,grid2,wt2,zeta1,zeta2,A12,B12,C12);
-      becke_weight_2d(gs,grid1,wtt1,grid2,wt2,Z1,Z2,A12,B12,C12);
+      becke_weight_2d(gs,grid1,wtt1,grid2,wt2,1,Z1,Z2,A12,B12,C12);
 
       //eliminate_small_wt(estart,gs,wtt1);
       //eliminate_small_wt(estart,gs,wt2);
diff --git a/src/integrals/integrals_aux.cpp b/src/integrals/integrals_aux.cpp
index d337609..5c566c1 100644
--- a/src/integrals/integrals_aux.cpp
+++ b/src/integrals/integrals_aux.cpp
@@ -1670,8 +1670,10 @@ vector<vector<double> > setup_integrals_gsgpu(vector<vector<double> >& basis_aux
 
 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 (S!=NULL)
+    get_overlap(S, N, natm, nbas, nenv, atm, bas, env);
+  if (jH1!=NULL)
+    get_hcore(jH1, N, natm, nbas, nenv, atm, bas, env);
   if (T!=NULL)
     get_tcore(T, N, natm, nbas, nenv, atm, bas, env);
 
@@ -1688,7 +1690,7 @@ void compute_integrals_g(int natm, int nbas, int nenv, int N, int Naux, int nbas
     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
+  else if (prl>0)
   {
     printf("  no auxiliary basis \n");
   }
@@ -1740,7 +1742,7 @@ void compute_integrals_g(int natm, int nbas, int nenv, int N, int Naux, int nbas
     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
+  else if (prl>0)
   {
     printf("  no auxiliary basis \n");
   }
diff --git a/src/libcintw/cintprep.cpp b/src/libcintw/cintprep.cpp
index d5c7f57..0b88ed9 100644
--- a/src/libcintw/cintprep.cpp
+++ b/src/libcintw/cintprep.cpp
@@ -15,7 +15,8 @@
 
 #define ROOTEIGHTPI 5.013256549262001
 
-CINTPrep::CINTPrep(bool doing_ri) {
+CINTPrep::CINTPrep(bool doing_ri)
+{
   atm = NULL;
   bas = NULL;
   env = NULL;
@@ -32,7 +33,8 @@ CINTPrep::CINTPrep(bool doing_ri) {
   nenv_ri = 0;
 }
 
-CINTPrep::~CINTPrep() {
+CINTPrep::~CINTPrep()
+{
   if (atm != NULL) {
     delete [] atm;
   }
@@ -49,37 +51,44 @@ double get_gto_norm(int shl1, double zeta1)
   return CINTgto_norm(shl1, zeta1)/ROOTEIGHTPI;
 }
 
-int CINTPrep::get_var_dim() {
+int CINTPrep::get_var_dim()
+{
   return var_dim;
 }
 
-int CINTPrep::get_nbas() {
+int CINTPrep::get_nbas()
+{
   return nbas;
 }
 
-int CINTPrep::get_nenv() {
+int CINTPrep::get_nenv()
+{
   return nenv;
 }
 
-int *CINTPrep::get_atm() {
+int *CINTPrep::get_atm()
+{
   int *tmp = atm;
   atm = NULL;
   return tmp;
 }
 
-int *CINTPrep::get_bas() {
+int *CINTPrep::get_bas()
+{
   int *tmp = bas;
   bas = NULL;
   return tmp;
 }
 
-double *CINTPrep::get_env() {
+double *CINTPrep::get_env()
+{
   double *tmp = env;
   env = NULL;
   return tmp;
 }
 
-void CINTPrep::set_atm(int *&atm_in) {
+void CINTPrep::set_atm(int *&atm_in)
+{
   if (atm != NULL) {
     delete [] atm;
   }
@@ -87,7 +96,8 @@ void CINTPrep::set_atm(int *&atm_in) {
   atm_in = NULL;
 }
 
-void CINTPrep::set_bas(int *&bas_in) {
+void CINTPrep::set_bas(int *&bas_in)
+{
   if (bas != NULL) {
     delete [] bas;
   }
@@ -95,7 +105,8 @@ void CINTPrep::set_bas(int *&bas_in) {
   bas_in = NULL;
 }
 
-void CINTPrep::set_env(double *&env_in) {
+void CINTPrep::set_env(double *&env_in)
+{
   if (env != NULL) {
     delete [] env;
   }
@@ -104,38 +115,42 @@ void CINTPrep::set_env(double *&env_in) {
 }
 
 //RI methods
-int CINTPrep::get_var_dim_ri() {
+int CINTPrep::get_var_dim_ri()
+{
   return var_dim_ri;
 }
 
-int CINTPrep::get_nbas_ri() {
+int CINTPrep::get_nbas_ri()
+{
   return nbas_ri;
 }
 
-int CINTPrep::get_nenv_ri() {
+int CINTPrep::get_nenv_ri()
+{
   return nenv_ri;
 }
 //end ri methods
 
-void CINTPrep::copy_atoms(vector< int > &atoms_copy) {
+void CINTPrep::copy_atoms(vector< int > &atoms_copy)
+{
   vector< int >().swap(atoms_copy);
   int size = atoms.size();
   atoms_copy.reserve(size);
-  for (int i = 0; i < size; i++) {
+  for (int i = 0; i < size; i++)
     atoms_copy.push_back(atoms[i]);
-  }
 }
 
-void CINTPrep::copy_coord(vector< double > &coord_copy) {
+void CINTPrep::copy_coord(vector< double > &coord_copy)
+{
   vector< double >().swap(coord_copy);
   int size = coord.size();
   coord_copy.reserve(size);
-  for (int i = 0; i < size; i++) {
+  for (int i = 0; i < size; i++)
     coord_copy.push_back(coord[i]);
-  }
 }
 
-void CINTPrep::assign_coords(int natoms, int *atomlist, double *coords, bool in_bohr) {
+void CINTPrep::assign_coords(int natoms, int *atomlist, double *coords, bool in_bohr)
+{
   vector< double >().swap(coord);
   vector< int >().swap(atoms);
   atoms.reserve(natoms);
@@ -143,15 +158,14 @@ void CINTPrep::assign_coords(int natoms, int *atomlist, double *coords, bool in_
 
   double scale = (in_bohr) ? 1. : ANG2BOHR;
 
-  for (int i = 0; i < natoms; i++) {
+  for (int i = 0; i < natoms; i++)
     atoms.push_back(atomlist[i]);
-  }
-  for (int i = 0; i < natoms * 3; i++) {
+  for (int i = 0; i < natoms * 3; i++)
     coord.push_back(coords[i]*scale);
-  }
 }
 
-void CINTPrep::assign_coords(int natoms, int *atomlist, float *coords, bool in_bohr) {
+void CINTPrep::assign_coords(int natoms, int *atomlist, float *coords, bool in_bohr)
+{
   vector< double >().swap(coord);
   vector< int >().swap(atoms);
   atoms.reserve(natoms);
@@ -159,24 +173,22 @@ void CINTPrep::assign_coords(int natoms, int *atomlist, float *coords, bool in_b
 
   double scale = (in_bohr) ? 1. : ANG2BOHR;
 
-  for (int i = 0; i < natoms; i++) {
+  for (int i = 0; i < natoms; i++)
     atoms.push_back(atomlist[i]);
-  }
-  for (int i = 0; i < natoms * 3; i++) {
+  for (int i = 0; i < natoms * 3; i++)
     coord.push_back(coords[i]*scale);
-  }
 }
 
-void CINTPrep::read_xyz(string inpxyz) {
+void CINTPrep::read_xyz(string inpxyz)
+{
   xyzfile = inpxyz;
   ifstream inpfile(xyzfile.c_str());
   vector< double >().swap(coord);
   vector< int >().swap(atoms);
   int natoms;
 
-  if (inpfile.fail()) {
+  if (inpfile.fail())
     printf("ERROR: could not open %s\n", xyzfile.c_str());
-  }
   printf("Reading coordinates from: %s\n", xyzfile.c_str());
 
   string line;
@@ -199,7 +211,8 @@ void CINTPrep::read_xyz(string inpxyz) {
   atoms.reserve(natoms);
   coord.reserve(natoms*3);
 
-  for (int i = 0; i < natoms; i++) {
+  for (int i = 0; i < natoms; i++)
+  {
     string symbol;
     int atomic_num;
     double x, y, z;
@@ -233,7 +246,8 @@ void CINTPrep::read_xyz(string inpxyz) {
   return;
 }
 
-void CINTPrep::read_bas(string inpbas) {
+void CINTPrep::read_bas(string inpbas)
+{
   if (basmap.size() > 0) {
     basmap.clear();
   }
@@ -254,12 +268,14 @@ void CINTPrep::read_bas(string inpbas) {
   string line;
   istringstream iss;
 
-  while (getline(inpfile, line)) {
+  while (getline(inpfile, line))
+  {
     lines.push_back(line);
   }
 
   int i = 0;
-  while (i < lines.size()) {
+  while (i < lines.size())
+  {
     if (!all_of(lines[i].begin(), lines[i].end(), [](char c) {return isspace(c);})
         && lines[i].c_str()[0] != '*') {
       iss.str(lines[i]);
@@ -274,11 +290,13 @@ void CINTPrep::read_bas(string inpbas) {
       i++;
 
       int atom_num = elem_2_int.at(atom_symbol);
-      if (present_atoms.count(atom_num)) {
+      if (present_atoms.count(atom_num))
+      {
         present_atoms.erase(atom_num);
         basis_t basis;
         basis.nuc = atom_num;
-        while (i < lines.size()) {
+        while (i < lines.size())
+        {
           if (!all_of(lines[i].begin(), lines[i].end(), [](char c) {return isspace(c);})
               && lines[i].c_str()[0] != '*') {
             iss.str(lines[i]);
@@ -379,7 +397,8 @@ void CINTPrep::read_bas(string inpbas) {
   }
 }
 
-bool CINTPrep::read_bas_ri(string auxbas) {
+bool CINTPrep::read_bas_ri(string auxbas)
+{
   if (basmap_ri.size() > 0) {
     basmap_ri.clear();
   }
@@ -403,9 +422,11 @@ bool CINTPrep::read_bas_ri(string auxbas) {
   }
 
   int i = 0;
-  while (i < lines.size()) {
-    if (!all_of(lines[i].begin(), lines[i].end(), [](char c) {return isspace(c);}) 
-        && lines[i].c_str()[0] != '*') {
+  while (i < lines.size())
+  {
+    if (!all_of(lines[i].begin(), lines[i].end(), [](char c) {return isspace(c);})
+        && lines[i].c_str()[0] != '*')
+    {
       iss.str(lines[i]);
       string atom_symbol;
       iss >> atom_symbol;
@@ -417,13 +438,16 @@ bool CINTPrep::read_bas_ri(string auxbas) {
       i++;
 
       int atom_num = elem_2_int.at(atom_symbol);
-      if (present_atoms.count(atom_num)) {
+      if (present_atoms.count(atom_num))
+      {
         present_atoms.erase(atom_num);
         basis_t basis;
         basis.nuc = atom_num;
-        while (i < lines.size()) {
-          if (!all_of(lines[i].begin(), lines[i].end(), [](char c) {return isspace(c);}) 
-              && lines[i].c_str()[0] != '*') {
+        while (i < lines.size())
+        {
+          if (!all_of(lines[i].begin(), lines[i].end(), [](char c) {return isspace(c);})
+              && lines[i].c_str()[0] != '*')
+          {
             iss.str(lines[i]);
             string ang;
             int n_prim;
@@ -500,10 +524,13 @@ bool CINTPrep::read_bas_ri(string auxbas) {
 
         basmap_ri[atom_num] = basis;
       } // if present_atoms.count
-      else {
-        while (i < lines.size()) {
+      else
+      {
+        while (i < lines.size())
+        {
           i++;
-          if (lines[i].c_str()[0] == '*') {
+          if (lines[i].c_str()[0] == '*')
+          {
             i++;
             break;
           } // if lines[i]
@@ -512,7 +539,8 @@ bool CINTPrep::read_bas_ri(string auxbas) {
 
       iss.clear();
     }
-    else {
+    else
+    {
       i++;
     }
   } // while i < lines.size (outer)
@@ -520,13 +548,12 @@ bool CINTPrep::read_bas_ri(string auxbas) {
   {
     printf("  did not find an auxiliary basis set \n");
     return 0;
-    //printf("ERROR: Missing basis set for some atoms\n");
-    //exit(1);
   }
   return 1;
 }
 
-void CINTPrep::prep_env() {
+void CINTPrep::prep_env()
+{
   if (atm != NULL) {
     delete [] atm;
   }
@@ -543,18 +570,23 @@ void CINTPrep::prep_env() {
   int offset = PTR_ENV_START;
   var_dim = 0;
 
-  for (int i = 0; i < natm; i++) {
+  for (int i = 0; i < natm; i++)
+  {
     int nshls = basmap.at(atoms[i]).shells.size();
-    for (int j = 0; j < nshls; j++) {
+    for (int j = 0; j < nshls; j++)
+    {
       nbas++;
     } // for j
   } // for i
 
   nbas_ri = 0;
-  if (do_ri) {
-    for (int i = 0; i < natm; i++) {
+  if (do_ri)
+  {
+    for (int i = 0; i < natm; i++)
+    {
       int nshls = basmap_ri.at(atoms[i]).shells.size();
-      for (int j = 0; j < nshls; j++) {
+      for (int j = 0; j < nshls; j++)
+      {
         nbas_ri++;
       } // for j
     } // for i
@@ -568,7 +600,8 @@ void CINTPrep::prep_env() {
   bas = new int[nbas_all * BAS_SLOTS]();
   atm = new int[natm * ATM_SLOTS]();
 
-  for (int i = 0; i < natm; i++) {
+  for (int i = 0; i < natm; i++)
+  {
     int atom_num = atoms[i];
     atm[i * ATM_SLOTS + 0] = atom_num;
     atm[i * ATM_SLOTS + 1] = offset;
@@ -578,11 +611,13 @@ void CINTPrep::prep_env() {
   } // for i
 
   int bas_num = 0;
-  for (int i = 0; i < natm; i++) {
+  for (int i = 0; i < natm; i++)
+  {
     int atom_num = atoms[i];
 
     int N_at = 0;
-    for (int j = 0; j < basmap.at(atom_num).shells.size(); j++) {
+    for (int j = 0; j < basmap.at(atom_num).shells.size(); j++)
+    {
       int nshls = basmap.at(atom_num).shells[j];
       int nprim = basmap.at(atom_num).exps[j].size();
       bas[ATOM_OF  +BAS_SLOTS*bas_num] = i;
@@ -602,10 +637,13 @@ void CINTPrep::prep_env() {
 
   } // for i
 
-  if (do_ri) {
-    for (int i = 0; i < natm; i++) {
+  if (do_ri)
+  {
+    for (int i = 0; i < natm; i++)
+    {
       int atom_num = atoms[i];
-      for (int j = 0; j < basmap_ri.at(atom_num).shells.size(); j++) {
+      for (int j = 0; j < basmap_ri.at(atom_num).shells.size(); j++)
+      {
         int nshls = basmap_ri.at(atom_num).shells[j];
         int nprim = basmap_ri.at(atom_num).exps[j].size();
         bas[ATOM_OF  +BAS_SLOTS*bas_num] = i;
@@ -631,38 +669,38 @@ void CINTPrep::prep_env() {
   nenv = offset;
   env = new double[nenv]();
   offset = PTR_ENV_START;
-  for (int i = 0; i < natm*3; i++) {
+  for (int i=0;i<natm*3;i++)
     env[offset++] = coord[i];
-  }
 
-  for (int i = 0; i < natm; i++) {
+  for (int i = 0; i < natm; i++)
+  {
     int atom_num = atoms[i];
-    for (int j = 0; j < basmap.at(atom_num).shells.size(); j++) {
+    for (int j = 0; j < basmap.at(atom_num).shells.size(); j++)
+    {
       int nexp = basmap.at(atom_num).exps[j].size();
       int ncof = basmap.at(atom_num).coef[j].size();
       int shll = basmap.at(atom_num).shells[j];
-      for (int k = 0; k < nexp; k++) {
+      for (int k = 0; k < nexp; k++)
         env[offset+k] = basmap.at(atom_num).exps[j][k];
-      } // for k
-      for (int k = 0; k < ncof; k++) {
+      for (int k = 0; k < ncof; k++)
         env[offset + nexp + k] = basmap.at(atom_num).coef[j][k] * CINTgto_norm(shll, env[offset+k]);
-      } // for k
       offset += nexp + ncof;
     } // for j
   } // for i
-  if (do_ri) {
-    for (int i = 0; i < natm; i++) {
+  if (do_ri)
+  {
+    for (int i = 0; i < natm; i++)
+    {
       int atom_num = atoms[i];
-      for (int j = 0; j < basmap_ri.at(atom_num).shells.size(); j++) {
+      for (int j = 0; j < basmap_ri.at(atom_num).shells.size(); j++)
+      {
         int nexp = basmap_ri.at(atom_num).exps[j].size();
         int ncof = basmap_ri.at(atom_num).coef[j].size();
         int shll = basmap_ri.at(atom_num).shells[j];
-        for (int k = 0; k < nexp; k++) {
+        for (int k = 0; k < nexp; k++)
           env[offset+k] = basmap_ri.at(atom_num).exps[j][k];
-        } // for k
-        for (int k = 0; k < ncof; k++) {
+        for (int k = 0; k < ncof; k++)
           env[offset + nexp + k] = basmap_ri.at(atom_num).coef[j][k] * CINTgto_norm(shll, env[offset+k]);
-        } // for k
         offset += nexp + ncof;
       } // for j
     } // for i
diff --git a/src/libcintw/cintwrapper.cpp b/src/libcintw/cintwrapper.cpp
index 1a5f565..5b1c9c3 100644
--- a/src/libcintw/cintwrapper.cpp
+++ b/src/libcintw/cintwrapper.cpp
@@ -24,11 +24,13 @@ using namespace std;
 
 bool BT::DO_CART = true;
 
-int binom(int n, int k) {
+int binom(int n, int k)
+{
   return (tgamma(n+1))/(tgamma(n+1-k)*tgamma(k+1));
 }
 
-int calc_di(int i, int *bas) {
+int calc_di(int i, int *bas)
+{
   int idx_i = 0;
   for (int j = 0; j < i; j++) {
     idx_i += CINTcgto_cart(j, bas);
@@ -268,7 +270,8 @@ void get_overlap_ri(double * overlap, int Naux,
 
 void get_hcore(double *hcore, double *En, double *T, int N,
                int natm, int nbas, int nenv,
-               int *atm, int *bas, double *env) {
+               int *atm, int *bas, double *env)
+{
   int idx_i = 0;
   int di, dj;
   int shls[2];
@@ -739,12 +742,12 @@ void gen_4c_overlap_m4(double* ovlp4, size_t N,
  #if USE_PARA_4C
   delete [] shell_offset;
  #endif
-  
+
   #else
     printf("\n\n\n ERROR: this SlaterGPU is not set up for 4c center integrals gbasis \n");
     exit(-1);
   #endif
-  
+
   return;
 }
 
@@ -871,36 +874,229 @@ void gen_4c_overlap(double* ovlp4, size_t N,
   return;
 }
 
-void gen_eri(double **eri, int N,
-             int natm, int nbas, int nenv,
-             int *atm, int *bas, double *env)
+void gen_eri_4b(double** eri, int N1, int N2, int N3, int N,
+             int natm, int nbas1, int nbas2, int nbas3, int nbas, int nenv,
+             int* atm, int* bas, double* env, int nomp)
 {
-  //int N2 = N*N;
+  int N4 = N-N1-N2-N3;
+  int nbas12 = nbas1+nbas2;
+  int nbas13 = nbas12+nbas3;
+  if (nomp<1) nomp = 1;
+  if (nomp>nbas3) nomp = nbas3;
+
+  int idxi = 0;
+  int idxj = 0;
+  //int idxk = 0;
+  CINTOpt* no_opt = NULL;
+
+ //buffer size
+  int mdi = 0; int mdj = 0; int mdk = 0; int mdl = 0;
+  for (int i=0;i<nbas1;i++)
+  {
+    int di = CINTcgto_cart(i,bas);
+    if (di>mdi) mdi = di;
+  }
+  for (int j=nbas1;j<nbas12;j++)
+  {
+    int dj = CINTcgto_cart(j,bas);
+    if (dj>mdj) mdj = dj;
+  }
+  for (int k=nbas12;k<nbas13;k++)
+  {
+    int dk = CINTcgto_cart(k,bas);
+    if (dk>mdk) mdk = dk;
+  }
+  for (int l=nbas13;l<nbas;l++)
+  {
+    int dl = CINTcgto_cart(l,bas);
+    if (dl>mdl) mdl = dl;
+  }
+
+  int sizeb = mdi*mdj*mdk*mdl;
+  double* bufall = new double[sizeb*nomp];
+  if (BT::DO_CART)
+  {
+    for (int i=0;i<nbas1;i++)
+    {
+      int di = CINTcgto_cart(i,bas);
+      for (int j=nbas1;j<nbas12;j++)
+      {
+        int dj = CINTcgto_cart(j,bas);
+       #pragma omp parallel for num_threads(nomp)
+        for (int k=nbas12;k<nbas13;k++)
+        {
+          //int tid = 0;
+          int tid = omp_get_thread_num();
+          double* buf = &bufall[sizeb*tid];
+
+          int idxk = 0;
+          for (int k1=nbas12;k1<k;k1++) idxk += CINTcgto_cart(k1,bas);
+          int idxl = 0;
+
+          int dk = CINTcgto_cart(k,bas);
+          for (int l=nbas13;l<nbas;l++)
+          {
+            int dl = CINTcgto_cart(l,bas);
+            int shls[4];
+            shls[0] = i;
+            shls[1] = j;
+            shls[2] = k;
+            shls[3] = l;
+            cint2e_cart(buf, shls, atm, natm, bas, nbas, env, no_opt);
+
+            for (int l1=0;l1<dl;l1++)
+            {
+              int ol = l1 + idxl;
+              for (int k1=0;k1<dk;k1++)
+              {
+                int ok = k1 + idxk;
+                for (int j1=0;j1<dj;j1++)
+                {
+                  int oj = j1 + idxj;
+                  for (int i1=0;i1<di;i1++)
+                  {
+                    int oi = i1 + idxi;
+                    int ij = oi * N2 + oj;
+                    int kl = ok * N4 + ol;
+                    int bdx = l1*(dk*dj*di) + k1*(dj*di) + j1*(di) + i1;
+                    // printf("oi %2d oj %2d ok %2d ol %2d ij %6d kl %6d\n",oi, oj, ok, ol,ij,kl);
+                    eri[ij][kl] = buf[bdx];
+                  } // for i1
+                } // for j1
+              } // for k1
+            } // for l1
+            idxl += dl;
+          } // for l
+          idxl = 0;
+          //idxk += dk;
+        } // for k
+        //idxk = 0;
+        idxj += dj;
+      } // for j
+      idxj = 0;
+      idxi += di;
+    } // for i
+  } // if BT::DO_CART
+  else
+  {
+    for (int i=0;i<nbas1;i++)
+    {
+      int di = CINTcgto_spheric(i,bas);
+      for (int j=nbas1;j<nbas12;j++)
+      {
+        int dj = CINTcgto_spheric(j,bas);
+       #pragma omp parallel for num_threads(nomp)
+        for (int k=nbas12;k<nbas13;k++)
+        {
+         //CPMZ need shell offsets
+          //int tid = 0;
+          int tid = omp_get_thread_num();
+          double* buf = &bufall[sizeb*tid];
+
+          int idxk = 0;
+          for (int k1=nbas12;k1<k;k1++) idxk += CINTcgto_spheric(k1,bas);
+
+          int idxl = 0;
+          int dk = CINTcgto_spheric(k,bas);
+          for (int l=nbas13;l<nbas;l++)
+          {
+            int dl = CINTcgto_spheric(l,bas);
+            int shls[4];
+            shls[0] = i;
+            shls[1] = j;
+            shls[2] = k;
+            shls[3] = l;
+            cint2e_sph(buf, shls, atm, natm, bas, nbas, env, no_opt);
+
+            for (int l1=0;l1<dl;l1++)
+            {
+              int ol = l1 + idxl;
+              for (int k1=0;k1<dk;k1++)
+              {
+                int ok = k1 + idxk;
+                for (int j1=0;j1<dj;j1++)
+                {
+                  int oj = j1 + idxj;
+                  for (int i1=0;i1<di;i1++)
+                  {
+                    int oi = i1 + idxi;
+                    int ij = oi * N2 + oj;
+                    int kl = ok * N4 + ol;
+                    int bdx = l1*(dk*dj*di) + k1*(dj*di) + j1*(di) + i1;
+                    // printf("oi %2d oj %2d ok %2d ol %2d ij %6d kl %6d\n",oi, oj, ok, ol,ij,kl);
+                    eri[ij][kl] = buf[bdx];
+                  } // for i1
+                } // for j1
+              } // for k1
+            } // for l1
+            idxl += dl;
+          } // for l
+          idxl = 0;
+          //idxk += dk;
+        } // for k
+        //idxk = 0;
+        idxj += dj;
+      } // for j
+      idxj = 0;
+      idxi += di;
+    } // for i
+  } // else
+
+  delete [] bufall;
+
+  return;
+}
+
+void gen_eri_2b(double** eri, int N1, int N,
+             int natm, int nbas1, int nbas, int nenv,
+             int* atm, int* bas, double* env, int nomp)
+{
+  int N2 = N-N1;
+  if (nomp<1) nomp = 1;
 
   int idxi = 0;
   int idxj = 0;
   int idxk = 0;
-  int idxl = 0;
-  int di = 0;
-  int dj = 0;
-  int dk = 0;
-  int dl = 0;
-  CINTOpt *no_opt = NULL;
+  CINTOpt* no_opt = NULL;
+
+ //buffer size
+  int mdi = 0; int mdk = 0;
+  for (int i=0;i<nbas1;i++)
+  {
+    int di = CINTcgto_cart(i,bas);
+    if (di>mdi) mdi = di;
+  }
+  for (int k=nbas1;k<nbas;k++)
+  {
+    int dk = CINTcgto_cart(k,bas);
+    if (dk>mdk) mdk = dk;
+  }
+
+  int sizeb = mdi*mdi*mdk*mdk;
+  double* bufall = new double[sizeb*nomp];
   if (BT::DO_CART)
   {
-    for (int i = 0; i < nbas; i++)
+    for (int i=0;i<nbas1;i++)
     {
-      di = CINTcgto_cart(i,bas);
-      for (int j = 0; j < nbas; j++)
+      int di = CINTcgto_cart(i,bas);
+      for (int j=0;j<nbas1;j++)
       {
-        dj = CINTcgto_cart(j,bas);
-        for (int k = 0; k < nbas; k++)
+        int dj = CINTcgto_cart(j,bas);
+       #pragma omp parallel for num_threads(nomp)
+        for (int k=nbas1;k<nbas;k++)
         {
-          dk = CINTcgto_cart(k,bas);
-          for (int l = 0; l < nbas; l++)
+          //int tid = 0;
+          int tid = omp_get_thread_num();
+          double* buf = &bufall[sizeb*tid];
+
+          int idxk = 0;
+          for (int k1=nbas1;k1<k;k1++) idxk += CINTcgto_cart(k1,bas);
+
+          int idxl = 0;
+          int dk = CINTcgto_cart(k,bas);
+          for (int l=nbas1;l<nbas;l++)
           {
-            dl = CINTcgto_cart(l,bas);
-            double * buf = new double[di*dj*dk*dl];
+            int dl = CINTcgto_cart(l,bas);
             int shls[4];
             shls[0] = i;
             shls[1] = j;
@@ -908,24 +1104,27 @@ void gen_eri(double **eri, int N,
             shls[3] = l;
             cint2e_cart(buf, shls, atm, natm, bas, nbas, env, no_opt);
 
-            for (int l1 = 0; l1 < dl; l1++) {
-              for (int k1 = 0; k1 < dk; k1++) {
-                for (int j1 = 0; j1 < dj; j1++) {
-                  for (int i1 = 0; i1 < di; i1++) {
+            for (int l1=0;l1<dl;l1++)
+            {
+              int ol = l1 + idxl;
+              for (int k1=0;k1<dk;k1++)
+              {
+                int ok = k1 + idxk;
+                for (int j1=0;j1<dj;j1++)
+                {
+                  int oj = j1 + idxj;
+                  for (int i1=0;i1<di;i1++)
+                  {
                     int oi = i1 + idxi;
-                    int oj = j1 + idxj;
-                    int ok = k1 + idxk;
-                    int ol = l1 + idxl;
-                    int ij = oi * N + oj;
-                    int kl = ok * N + ol;
+                    int ij = oi * N1 + oj;
+                    int kl = ok * N2 + ol;
                     int bdx = l1*(dk*dj*di) + k1*(dj*di) + j1*(di) + i1;
                     // printf("oi %2d oj %2d ok %2d ol %2d ij %6d kl %6d\n",oi, oj, ok, ol,ij,kl);
-                    eri[ij][kl] = eri[kl][ij] = buf[bdx];
+                    eri[ij][kl] = buf[bdx];
                   } // for i1
                 } // for j1
               } // for k1
             } // for l1
-            delete [] buf;
             idxl += dl;
           } // for l
           idxl = 0;
@@ -940,19 +1139,27 @@ void gen_eri(double **eri, int N,
   } // if BT::DO_CART
   else
   {
-    for (int i = 0; i < nbas; i++)
+    for (int i=0;i<nbas1;i++)
     {
-      di = CINTcgto_spheric(i,bas);
-      for (int j = 0; j < nbas; j++)
+      int di = CINTcgto_spheric(i,bas);
+      for (int j=0;j<nbas1;j++)
       {
-        dj = CINTcgto_spheric(j,bas);
-        for (int k = 0; k < nbas; k++)
+        int dj = CINTcgto_spheric(j,bas);
+       #pragma omp parallel for num_threads(nomp)
+        for (int k=nbas1;k<nbas;k++)
         {
-          dk = CINTcgto_spheric(k,bas);
-          for (int l = 0; l < nbas; l++)
+          //int tid = 0;
+          int tid = omp_get_thread_num();
+          double* buf = &bufall[sizeb*tid];
+
+          int idxk = 0;
+          for (int k1=nbas1;k1<k;k1++) idxk += CINTcgto_spheric(k1,bas);
+
+          int idxl = 0;
+          int dk = CINTcgto_spheric(k,bas);
+          for (int l=nbas1;l<nbas;l++)
           {
-            dl = CINTcgto_spheric(l,bas);
-            double * buf = new double[di*dj*dk*dl];
+            int dl = CINTcgto_spheric(l,bas);
             int shls[4];
             shls[0] = i;
             shls[1] = j;
@@ -960,24 +1167,27 @@ void gen_eri(double **eri, int N,
             shls[3] = l;
             cint2e_sph(buf, shls, atm, natm, bas, nbas, env, no_opt);
 
-            for (int l1 = 0; l1 < dl; l1++) {
-              for (int k1 = 0; k1 < dk; k1++) {
-                for (int j1 = 0; j1 < dj; j1++) {
-                  for (int i1 = 0; i1 < di; i1++) {
+            for (int l1=0;l1<dl;l1++)
+            {
+              int ol = l1 + idxl;
+              for (int k1=0;k1<dk;k1++)
+              {
+                int ok = k1 + idxk;
+                for (int j1=0;j1<dj;j1++)
+                {
+                  int oj = j1 + idxj;
+                  for (int i1=0;i1<di;i1++)
+                  {
                     int oi = i1 + idxi;
-                    int oj = j1 + idxj;
-                    int ok = k1 + idxk;
-                    int ol = l1 + idxl;
-                    int ij = oi * N + oj;
-                    int kl = ok * N + ol;
+                    int ij = oi * N1 + oj;
+                    int kl = ok * N2 + ol;
                     int bdx = l1*(dk*dj*di) + k1*(dj*di) + j1*(di) + i1;
                     // printf("oi %2d oj %2d ok %2d ol %2d ij %6d kl %6d\n",oi, oj, ok, ol,ij,kl);
-                    eri[ij][kl] = eri[kl][ij] = buf[bdx];
+                    eri[ij][kl] = buf[bdx];
                   } // for i1
                 } // for j1
               } // for k1
             } // for l1
-            delete [] buf;
             idxl += dl;
           } // for l
           idxl = 0;
@@ -991,12 +1201,245 @@ void gen_eri(double **eri, int N,
     } // for i
   } // else
 
+  delete [] bufall;
+
+  return;
+}
+
+void gen_eri_2(double** eri, int N,
+             int natm, int nbas, int nenv,
+             int* atm, int* bas, double* env, int nomp)
+{
+  CINTOpt* no_opt = NULL;
+  if (nomp<1) nomp = 1;
+
+  if (BT::DO_CART)
+  { printf("\n ERROR: gen_eri_2 does not support Cartesian GTOs \n"); exit(-1); }
+
+  int mdi = 0;
+  for (int i=0;i<nbas;i++)
+  {
+    int di = CINTcgto_cart(i,bas);
+    if (di>mdi) mdi = di;
+  }
+
+  int sbuf = mdi*mdi*mdi*mdi;
+  double* bufall = new double[sbuf*nomp];
+
+  vector<int> nf(nbas);
+  vector<int> offs(nbas);
+  int cum = 0;
+  for (int i=0;i<nbas;i++)
+  {
+    nf[i] = CINTcgto_spheric(i, bas);
+    offs[i] = cum; cum += nf[i];
+  }
+
+  //unique shell pairs
+  struct SP { int i, j; };
+  vector<SP> pairs;
+  pairs.reserve(nbas*(nbas+1)/2);
+  for (int i=0;i<nbas;i++)
+  for (int j=0;j<=i;j++)
+    pairs.push_back({i, j});
+  int npair = (int)pairs.size();
+
+  if (nomp>npair) nomp = npair;
+
+  int shls[4];
+  #pragma omp parallel for num_threads(nomp) schedule(dynamic) private(shls)
+  for (int p=0;p<npair;p++)
+  {
+    int tid = omp_get_thread_num();
+    double* buf = &bufall[tid*sbuf];
+
+    int i = pairs[p].i, j = pairs[p].j;
+    int di = nf[i], dj = nf[j], idxi = offs[i], idxj = offs[j];
+
+    for (int q=0;q<=p;q++) // ket pair index <= bra pair index -> bra-ket swap folded in
+    {
+      int k = pairs[q].i; int l = pairs[q].j;
+      int dk = nf[k]; int dl = nf[l]; int idxk = offs[k]; int idxl = offs[l];
+
+      shls[0] = i; shls[1] = j; shls[2] = k; shls[3] = l;
+      cint2e_sph(buf,shls,atm,natm,bas,nbas,env,no_opt);
+
+      for (int l1=0;l1<dl;l1++)
+      for (int k1=0;k1<dk;k1++)
+      for (int j1=0;j1<dj;j1++)
+      for (int i1=0;i1<di;i1++)
+      {
+        int oi = i1+idxi; int oj = j1+idxj; int ok = k1+idxk; int ol = l1+idxl;
+        int bdx = l1*(dk*dj*di) + k1*(dj*di) + j1*di + i1;
+        double val = buf[bdx];
+
+        int ij = oi*N+oj; int ji = oj*N+oi;
+        int kl = ok*N+ol; int lk = ol*N+ok;
+        eri[ij][kl]=val; eri[ji][kl]=val; eri[ij][lk]=val; eri[ji][lk]=val;
+        eri[kl][ij]=val; eri[kl][ji]=val; eri[lk][ij]=val; eri[lk][ji]=val;
+      }
+    }
+  }
+
+  delete [] bufall;
+  return;
+}
+
+void gen_eri(double** eri, int N,
+             int natm, int nbas, int nenv,
+             int* atm, int* bas, double* env, int nomp)
+{
+  CINTOpt* no_opt = NULL;
+  if (nomp<1) nomp = 1;
+
+  int mdi = 0;
+  for (int i=0;i<nbas;i++)
+  {
+    int di = CINTcgto_cart(i,bas);
+    if (di>mdi) mdi = di;
+  }
+
+  int sbuf = mdi*mdi*mdi*mdi;
+  double* bufall = new double[sbuf*nomp];
+
+  if (BT::DO_CART)
+  {
+    int shls[4];
+    int idxi = 0;
+    for (int i=0;i<nbas;i++)
+    {
+      int di = CINTcgto_cart(i,bas);
+
+     #pragma omp parallel for num_threads(nomp) private(shls)
+      for (int j=0;j<nbas;j++)
+      {
+        int tid = omp_get_thread_num();
+        double* buf = &bufall[tid*sbuf];
+
+        int idxj = 0; for (int j1=0;j1<j;j1++) idxj += CINTcgto_cart(j1,bas);
+
+        int dj = CINTcgto_cart(j,bas);
+        int idxk = 0;
+        for (int k=0;k<nbas;k++)
+        {
+          int dk = CINTcgto_cart(k,bas);
+          int idxl = 0;
+          for (int l=0;l<nbas;l++)
+          {
+            int dl = CINTcgto_cart(l,bas);
+            shls[0] = i;
+            shls[1] = j;
+            shls[2] = k;
+            shls[3] = l;
+            cint2e_cart(buf, shls, atm, natm, bas, nbas, env, no_opt);
+
+            for (int l1=0;l1<dl;l1++)
+            for (int k1=0;k1<dk;k1++)
+            for (int j1=0;j1<dj;j1++)
+            for (int i1=0;i1<di;i1++)
+            {
+              int oi = i1 + idxi;
+              int oj = j1 + idxj;
+              int ok = k1 + idxk;
+              int ol = l1 + idxl;
+              int ij = oi * N + oj;
+              int kl = ok * N + ol;
+              int bdx = l1*(dk*dj*di) + k1*(dj*di) + j1*(di) + i1;
+              // printf("oi %2d oj %2d ok %2d ol %2d ij %6d kl %6d\n",oi, oj, ok, ol,ij,kl);
+              eri[ij][kl] = buf[bdx];
+              //eri[ij][kl] = eri[kl][ij] = buf[bdx];
+
+            } // for l1k1j1l1
+            idxl += dl;
+          } // for l
+          idxl = 0;
+          idxk += dk;
+        } // for k
+        idxk = 0;
+        //idxj += dj;
+      } // for j
+      //idxj = 0;
+      idxi += di;
+    } // for i
+  } // if BT::DO_CART
+  else
+  {
+    int shls[4];
+    int idxi = 0;
+    for (int i=0;i<nbas;i++)
+    {
+      int di = CINTcgto_spheric(i,bas);
+     #pragma omp parallel for num_threads(nomp) private(shls)
+      for (int j=0;j<nbas;j++)
+      {
+        int tid = omp_get_thread_num();
+        double* buf = &bufall[tid*sbuf];
+
+        int idxj = 0; for (int j1=0;j1<j;j1++) idxj += CINTcgto_spheric(j1,bas);
+
+        int dj = CINTcgto_spheric(j,bas);
+
+        int idxk = 0;
+        for (int k=0;k<nbas;k++)
+        {
+          int dk = CINTcgto_spheric(k,bas);
+
+          int idxl = 0;
+          for (int l=0;l<nbas;l++)
+          {
+            int dl = CINTcgto_spheric(l,bas);
+            shls[0] = i;
+            shls[1] = j;
+            shls[2] = k;
+            shls[3] = l;
+            cint2e_sph(buf, shls, atm, natm, bas, nbas, env, no_opt);
+
+            for (int l1=0;l1<dl;l1++)
+            for (int k1=0;k1<dk;k1++)
+            for (int j1=0;j1<dj;j1++)
+            for (int i1=0;i1<di;i1++)
+            {
+              int oi = i1 + idxi;
+              int oj = j1 + idxj;
+              int ok = k1 + idxk;
+              int ol = l1 + idxl;
+              int ij = oi * N + oj;
+              int kl = ok * N + ol;
+              int bdx = l1*(dk*dj*di) + k1*(dj*di) + j1*(di) + i1;
+              // printf("oi %2d oj %2d ok %2d ol %2d ij %6d kl %6d\n",oi, oj, ok, ol,ij,kl);
+              eri[ij][kl] = buf[bdx];
+              //eri[ij][kl] = eri[kl][ij] = buf[bdx];
+
+            } // for l1k1j1i1
+            idxl += dl;
+          } // for l
+          idxl = 0;
+          idxk += dk;
+        } // for k
+        idxk = 0;
+        //idxj += dj;
+      } // for j
+      //idxj = 0;
+      idxi += di;
+    } // for i
+  } // else
+
+  delete [] bufall;
+
   return;
 }
 
-void gen_jMOI_gto(double **eri, int N,
+void gen_eri(double** eri, int N,
+             int natm, int nbas, int nenv,
+             int* atm, int* bas, double* env)
+{
+  return gen_eri(eri,N,natm,nbas,nenv,atm,bas,env,1);
+}
+
+void gen_jMOI_gto(double** eri, int N,
                   int natm, int nbas, int nenv,
-                  int *atm, int *bas, double *env) {
+                  int *atm, int *bas, double *env)
+{
   int N2 = N*N;
 
   int idxi = 0;
@@ -1008,7 +1451,8 @@ void gen_jMOI_gto(double **eri, int N,
   int dk = 0;
   int dl = 0;
   CINTOpt *no_opt = NULL;
-  if (BT::DO_CART) {
+  if (BT::DO_CART)
+  {
     for (int i = 0; i < nbas; i++) {
       di = CINTcgto_cart(i,bas);
       for (int j = 0; j < nbas; j++) {
@@ -1057,7 +1501,8 @@ void gen_jMOI_gto(double **eri, int N,
       idxi += di;
     } // for i
   } // if BT::DO_CART
-  else {
+  else
+  {
     for (int i = 0; i < nbas; i++) {
       di = CINTcgto_spheric(i,bas);
       for (int j = 0; j < nbas; j++) {
@@ -1301,17 +1746,21 @@ void gen_eri_3c(double *eri, int N, int Naux,
 
 void gen_eri_2c(double *eri, int Naux,
                 int natm, int nbas, int nenv, int nbas_ri,
-                int *atm, int *bas, double *env) {
+                int *atm, int *bas, double *env)
+{
   int Naux2 = Naux * Naux;
   int idxmu = 0;
   int idxnu = 0;
   int dmu = 0;
   int dnu = 0;
   CINTOpt *no_opt = NULL;
-  if (BT::DO_CART) {
-    for (int mu = 0; mu < nbas_ri; mu++) {
+  if (BT::DO_CART)
+  {
+    for (int mu = 0; mu < nbas_ri; mu++)
+    {
       dmu = CINTcgto_cart(mu+nbas,bas);
-      for (int nu = 0; nu <= mu; nu++) {
+      for (int nu = 0; nu <= mu; nu++)
+      {
         dnu = CINTcgto_cart(nu+nbas,bas);
         double *buf = new double[dmu * dnu];
         int shls[4];
@@ -1320,8 +1769,10 @@ void gen_eri_2c(double *eri, int Naux,
         shls[2] = nu + nbas;
         shls[3] = nbas + nbas_ri;
         cint2e_cart(buf, shls, atm, natm, bas, nbas + nbas_ri + 1, env, no_opt);
-        for (int mu1 = 0; mu1 < dmu; mu1++) {
-          for (int nu1 = 0; nu1 < dnu; nu1++) {
+        for (int mu1 = 0; mu1 < dmu; mu1++)
+        {
+          for (int nu1 = 0; nu1 < dnu; nu1++)
+          {
             int omu = mu1 + idxmu;
             int onu = nu1 + idxnu;
             int bdx = nu1 * dmu + mu1;
diff --git a/src/libio/CMakeLists.txt b/src/libio/CMakeLists.txt
index 31b3aae..1dd7d2e 100644
--- a/src/libio/CMakeLists.txt
+++ b/src/libio/CMakeLists.txt
@@ -2,6 +2,7 @@ set(IO_SOURCES
     read.cpp
     write.cpp
     print.cpp
+    molden.cpp
     )
 
 add_library(io SHARED "${IO_SOURCES}")
diff --git a/src/libio/molden.cpp b/src/libio/molden.cpp
new file mode 100644
index 0000000..dbaa572
--- /dev/null
+++ b/src/libio/molden.cpp
@@ -0,0 +1,736 @@
+#include "write.h"
+#include <sstream>
+
+bool close_val(double v1, double v2);
+
+void get_nxyzr(int n1, int l1, int m1, int& nx, int& ny, int& nz, int& nr)
+{
+  nx=ny=nz=nr=0;
+  if (n1==1)
+  {
+  }
+  else if (n1==2)
+  {
+    if (l1==0)
+      nr = 1;
+    else if (l1==1)
+    {
+      if (m1==0)      nz = 1;
+      else if (m1==1) nx = 1;
+      else            ny = 1;
+    }
+  }
+  else if (n1==3)
+  {
+    if (l1==0)
+      nr = 2;
+    else if (l1==1)
+    {
+      nr = 1;
+      if (m1==0)      nz = 1;
+      else if (m1==1) nx = 1;
+      else            ny = 1;
+    }
+    else if (l1==2)
+    {
+      nr = 0;
+     #if CART_D
+      if (m1==0) { nx = 2; }
+      if (m1==1) { ny = 2; }
+      if (m1==2) { nz = 2; }
+      if (m1==3) { nx = 1; ny = 1; }
+      if (m1==4) { nx = 1; nz = 1; }
+      if (m1==5) { ny = 1; nz = 1; }
+     #else
+      if (m1==-2) { nx = 1; ny = 1; }
+      if (m1==-1) { ny = 1; nz = 1; }
+      if (m1== 0) { nz = 2; }
+      if (m1== 1) { nx = 1; nz = 1; }
+      if (m1== 2) { nx = 2; ny = 2; }
+     #endif
+    }
+  }
+  else if (n1==4)
+  {
+    //printf(" WARNING: n=4 MO print incomplete \n");
+    if (l1==0)
+      nr = 3;
+    else if (l1==1)
+    {
+      nr = 2;
+      if (m1==0)      nz = 1;
+      else if (m1==1) nx = 1;
+      else            ny = 1;
+    }
+    else if (l1==2)
+    {
+      nr = 1;
+     #if CART_D
+      if (m1==0) { nx = 2; }
+      if (m1==1) { ny = 2; }
+      if (m1==2) { nz = 2; }
+      if (m1==3) { nx = 1; ny = 1; }
+      if (m1==4) { nx = 1; nz = 1; }
+      if (m1==5) { ny = 1; nz = 1; }
+     #else
+      if (m1==-2) { nx = 1; ny = 1; }
+      if (m1==-1) { ny = 1; nz = 1; }
+      if (m1== 0) { nz = 2; }
+      if (m1== 1) { nx = 1; nz = 1; }
+      if (m1== 2) { nx = 2; ny = 2; }
+     #endif
+    }
+    else if (l1==3)
+    {
+     //INCOMPLETE
+      if (m1==-2) { nx = 1; ny = 1; nz = 1; }
+      else if (m1==0) { nz = 3; }
+      else
+        nx = nx = nz = 444;
+    }
+
+  }
+  return;
+}
+
+void write_molden_mo_g(ofstream& outfile, vector<vector<double> > &basis, double* jCA, int No, double* eig)
+{
+  int N = basis.size();
+
+  outfile << "[MO]" << endl;
+  for (int i=0;i<N;i++)
+  {
+    int occ = 0; if (i<No) occ = 1;
+
+    outfile << "Sym=X" << endl;
+    if (eig!=NULL)
+      outfile << "Ene=" << eig[i] << endl;
+    else
+      outfile << "Ene=0." << endl;
+    outfile << "Spin=Alpha" << endl;
+    outfile << "Occup=" << occ << endl;
+
+    for (int j=0;j<N;j++)
+    {
+      int l1 = basis[j][1];
+      if (l1<2)
+        outfile << " " << j+1 << "  " << jCA[j*N+i] << endl;
+      else if (l1==2) //d orbs
+      {
+        int j1 = j+2; //0
+        outfile << " " << j+1 << "  " << jCA[j1*N+i] << endl;
+        j1 = j+3; //+1
+        outfile << " " << j+2 << "  " << jCA[j1*N+i] << endl;
+        j1 = j+1; //-1
+        outfile << " " << j+3 << "  " << jCA[j1*N+i] << endl;
+        j1 = j+4; //+2
+        outfile << " " << j+4 << "  " << jCA[j1*N+i] << endl;
+        j1 = j; //-2
+        outfile << " " << j+5 << "  " << jCA[j1*N+i] << endl;
+
+        j += 4; //jump ahead
+      }
+      else if (l1==3) //f orbs
+      {
+        int j1 = j+3; //0
+        outfile << " " << j+1 << "  " << jCA[j1*N+i] << endl;
+        j1 = j+4; //1
+        outfile << " " << j+2 << "  " << jCA[j1*N+i] << endl;
+        j1 = j+2; //-1
+        outfile << " " << j+3 << "  " << jCA[j1*N+i] << endl;
+        j1 = j+5; //2
+        outfile << " " << j+4 << "  " << jCA[j1*N+i] << endl;
+        j1 = j+1; //-2
+        outfile << " " << j+5 << "  " << jCA[j1*N+i] << endl;
+        j1 = j+6; //3
+        outfile << " " << j+6 << "  " << jCA[j1*N+i] << endl;
+        j1 = j; //-3
+        outfile << " " << j+7 << "  " << jCA[j1*N+i] << endl;
+
+        j += 6;
+      }
+      else if (l1==4) //g orbs
+      {
+        int j1 = j+4;
+        outfile << " " << j+1 << "  " << jCA[j1*N+i] << endl;
+        j1 = j+5;
+        outfile << " " << j+2 << "  " << jCA[j1*N+i] << endl;
+        j1 = j+3;
+        outfile << " " << j+3 << "  " << jCA[j1*N+i] << endl;
+        j1 = j+6;
+        outfile << " " << j+4 << "  " << jCA[j1*N+i] << endl;
+        j1 = j+2;
+        outfile << " " << j+5 << "  " << jCA[j1*N+i] << endl;
+        j1 = j+7;
+        outfile << " " << j+6 << "  " << jCA[j1*N+i] << endl;
+        j1 = j+1;
+        outfile << " " << j+7 << "  " << jCA[j1*N+i] << endl;
+        j1 = j+8;
+        outfile << " " << j+8 << "  " << jCA[j1*N+i] << endl;
+        j1 = j;
+        outfile << " " << j+9 << "  " << jCA[j1*N+i] << endl;
+
+        j += 8;
+      }
+      else //not printing other orbs at the moment
+        outfile << " " << j+1 << "  0.0" << endl;
+    }
+  }
+
+  int lmax = 0;
+  for (int i=0;i<N;i++)
+    lmax = max(lmax,(int)basis[i][1]);
+
+  if (lmax==2)
+    outfile << " [5D] " << endl;
+  else if (lmax==3)
+    outfile << " [5D7F] " << endl;
+  else
+    outfile << " [5D7F9G] " << endl;
+
+  return;
+}
+
+void write_molden_g(int btype, int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, double* eig, string fname)
+{
+ //this ftn expects orbitals to be spherical, not Cartesian. May crash for Cartesians
+  int N = basis.size();
+
+  No = abs(No);
+
+  const double n0 = norm_sh(0,0);
+  double B2A = 1./A2B;
+
+  string filename = fname;
+  ofstream outfile;
+  outfile.open(filename.c_str());
+  outfile << fixed << setprecision(8);
+
+  outfile << "[Molden Format]" << endl;
+  outfile << "[Atoms] (Angs)" << endl;
+  for (int i=0;i<natoms;i++)
+  {
+    string aname = get_aname(atno[i]);
+    double x1 = coords[3*i]*B2A; double y1 = coords[3*i+1]*B2A; double z1 = coords[3*i+2]*B2A;
+    outfile << " " << aname << "     " << i+1 << "   " << atno[i] << "   ";
+    outfile << x1 << "   " << y1 << "   " << z1 << endl;
+  }
+
+
+  outfile << "[GTO]" << endl;
+
+ //get molden formatted basis from basis vector
+  for (int i=0;i<natoms;i++)
+  {
+    string aname = get_aname(atno[i]);
+
+    int wb1 = -1;
+    for (int j=0;j<basis.size();j++)
+    if (basis[j][9]==i)
+    { wb1 = j; break; }
+    int wb2 = wb1;
+    for (int j=wb1+1;j<basis.size();j++)
+    if (basis[j][9]!=i)
+    { wb2 = j; break; }
+    if (i+1==natoms) wb2 = N;
+    //printf(" wb12: %2i %2i \n",wb1,wb2);
+
+    string b1 = "";
+    for (int j=wb1;j<wb2;)
+    {
+      int l1 = (int)basis[j][1];
+      int ng = (int)basis[j][3];
+      int row0 = j + l1; // the m==0 row of this shell
+      double nsh0 = norm_sh(l1,0);
+
+      string lv = "S"; if (l1==1) lv = "P"; else if (l1==2) lv = "D"; else if (l1==3) lv = "F"; else if (l1==4) lv = "G"; else if (l1==5) lv = "H";
+      b1 += lv+" "+to_string(ng)+" 1.0 \n";
+      for (int k=0;k<ng;k++)
+      {
+        double g1 = basis[row0][10+k];
+        double n1 = basis[row0][10+ng+k] * n0 / ( get_gto_norm(l1,g1) * nsh0 );
+        b1 += " " + to_string(g1) + "  " + to_string(n1) + "\n";
+      }
+
+      j += (2*l1+1);
+    }
+    outfile << i+1 << "    0 " << endl;
+    outfile << b1 << endl;
+  }
+
+  write_molden_mo_g(outfile,basis,jCA,No,eig);
+
+  outfile.close();
+
+  return;
+}
+
+void write_molden_g(int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, double* eig, string fname)
+{
+  return write_molden_g(0,natoms,atno,coords,basis,jCA,No,eig,fname);
+}
+
+void write_molden_ss(int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, string fname)
+{
+  int N = basis.size();
+
+  string filename = fname;
+  ofstream outfile;
+  outfile.open(filename.c_str());
+
+  double B2A = 1./A2B;
+
+  outfile << fixed << setprecision(8);
+
+  outfile << "[Molden Format]" << endl;
+  outfile << "[Atoms] (Angs)" << endl;
+  for (int i=0;i<natoms;i++)
+  {
+    string aname = get_aname(atno[i]);
+    double x1 = coords[3*i]*B2A; double y1 = coords[3*i+1]*B2A; double z1 = coords[3*i+2]*B2A;
+    outfile << " " << aname << "     " << i+1 << "   " << atno[i] << "   ";
+    outfile << x1 << "   " << y1 << "   " << z1 << endl;
+  }
+
+
+  outfile << "[STO]" << endl;
+
+  int nss = 4;
+  for (int i=0;i<natoms;i++)
+  {
+    double x1 = coords[3*i]; double y1 = coords[3*i+1]; double z1 = coords[3*i+2];
+
+    for (int j=0;j<N;j++)
+    {
+      if (close_val(x1,basis[j][5]) && close_val(y1,basis[j][6]) && close_val(z1,basis[j][7]))
+      {
+        int n1 = basis[j][0]; int l1 = basis[j][1]; int m1 = basis[j][2]; double zeta = basis[j][3]*B2A; double norm = basis[j][4];
+        int nx,ny,nz,nr;
+        get_nxyzr(n1,l1,m1,nx,ny,nz,nr);
+
+        if (l1<2)
+        {
+          for (int n=0;n<nss;n++)
+          {
+           //nr+n-l1?
+            double norm1 = norm*pow(zeta,n); if (n==2) norm1 *= 0.4; if (n==3) norm1 /= 15.;
+            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+            outfile << nr+n << " " << zeta << " " << norm1 << " " << endl;
+          }
+        }
+        else if (l1==2)
+        {
+         //for 2 of the 5 3d functions, divide into parts
+          if (nz==2)
+          {
+           //3dz2 --> 2.z2 - x2 - y2
+            for (int n=0;n<nss;n++)
+            {
+              double norm1 = norm*pow(zeta,n); if (n==2) norm1 *= 0.4; if (n==3) norm1 /= 15.;
+              nz = 2; nx = ny = 0;
+              outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+              outfile << nr+n << " " << zeta << " " << 2.*norm1 << " " << endl;
+              nz = 0; nx = 2; ny = 0;
+              outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+              outfile << nr+n << " " << zeta << " " << norm1 << " " << endl;
+              nz = 0; nx = 0; ny = 2;
+              outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+              outfile << nr+n << " " << zeta << " " << norm1 << " " << endl;
+            }
+          }
+          else if (nx==2 && ny==2)
+          {
+            //x2-y2 --> x2 - y2
+            for (int n=0;n<nss;n++)
+            {
+              double norm1 = norm*pow(zeta,n); if (n==2) norm1 *= 0.4; if (n==3) norm1 /= 15.;
+              nx = 2; ny = 0;
+              outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+              outfile << nr+n << " " << zeta << " " << norm1 << " " << endl;
+              nx = 0; ny = 2;
+              outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+              outfile << nr+n << " " << zeta << " " << norm1 << " " << endl;
+            }
+          }
+          else
+          {
+            for (int n=0;n<nss;n++)
+            {
+              double norm1 = norm*pow(zeta,n); if (n==2) norm1 *= 0.4; if (n==3) norm1 /= 15.;
+              outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+              outfile << nr+n << " " << zeta << " " << norm1 << " " << endl;
+            }
+          }
+        } //if l1==2
+
+       #if 0
+       //just drop most f or higher functions
+        else if (l1==3)
+        {
+         //m==-2
+          if (nx==1 && ny==1 && nz==1) //fxyz is the simplest 4f orbital
+          {
+            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+            outfile << nr << " " << zeta << " " << norm << " " << endl;
+          }
+         //m==0
+          if (nz==3)
+          {
+           //5z3
+            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+            outfile << nr << " " << zeta << " " << 5.*norm << " " << endl;
+           //-3zr2
+            nz = 1; nr = 2;
+            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+            outfile << nr << " " << zeta << " " << -3.*norm << " " << endl;
+          }
+          else
+          {
+            //incomplete
+          }
+        } //if l1==3
+       #endif
+
+      }
+    }
+  }
+
+  outfile << "[MO]" << endl;
+  for (int i=0;i<N;i++)
+  {
+    int occ = 0; if (i<No) occ = 1;
+
+    outfile << "Sym=X" << endl;
+    outfile << "Ene=0.0" << endl;
+    outfile << "Spin=Alpha" << endl;
+    outfile << "Occup=" << occ << endl;
+
+    int wj = 0;
+    for (int j=0;j<N;j++)
+    {
+      int n1 = basis[j][0]; int l1 = basis[j][1]; int m1 = basis[j][2];
+      int nx,ny,nz,nr;
+      get_nxyzr(n1,l1,m1,nx,ny,nz,nr);
+
+      if (l1<2)
+      {
+        for (int n=0;n<nss;n++)
+          outfile << " " << 1+wj++ << "  " << jCA[j*N+i] << endl;
+      }
+      else if (l1==2)
+      {
+        if (nz==2)
+        {
+          for (int n=0;n<nss;n++)
+          {
+            outfile << " " << 1+wj++ << "  " <<  jCA[j*N+i] << endl;
+            outfile << " " << 1+wj++ << "  " << -jCA[j*N+i] << endl;
+            outfile << " " << 1+wj++ << "  " << -jCA[j*N+i] << endl;
+          }
+        }
+        else if (nx==2 && ny==2)
+        {
+          for (int n=0;n<nss;n++)
+          {
+            outfile << " " << 1+wj++ << "  " <<  jCA[j*N+i] << endl;
+            outfile << " " << 1+wj++ << "  " << -jCA[j*N+i] << endl;
+          }
+        }
+        else
+        {
+          for (int n=0;n<nss;n++)
+            outfile << " " << 1+wj++ << "  " << jCA[j*N+i] << endl;
+        }
+      }
+     #if 0
+      else if (l1==3)
+      {
+       //incomplete
+        if (nx==1 && ny==1 && nz==1) //m==-2
+          outfile << " " << 1+wj++ << "  " << jCA[j*N+i] << endl;
+        else if (nz==3) //m==0
+          outfile << " " << 1+wj++ << "  " << jCA[j*N+i] << endl;
+      }
+     #endif
+    }
+  }
+
+  outfile.close();
+
+  return;
+}
+
+void write_molden(bool gbasis, int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, double* eig, string fname)
+{
+  if (gbasis) return write_molden_g(natoms,atno,coords,basis,jCA,No,eig,fname);
+  if (basis[0].size()>10) return write_molden_ss(natoms,atno,coords,basis,jCA,No,fname);
+
+  int N = basis.size();
+
+  bool missing_ftns = 0;
+  for (int j=0;j<N;j++)
+  if (basis[j][0]>3)
+    missing_ftns = 1;
+  if (missing_ftns)
+    printf(" WARNING: n=4 MO print incomplete \n");
+  //printf("\n WARNING: no f function printing to molden \n");
+
+  string filename = fname;
+  ofstream outfile;
+  outfile.open(filename.c_str());
+
+  double B2A = 1./A2B;
+
+  outfile << fixed << setprecision(8);
+
+  outfile << "[Molden Format]" << endl;
+  outfile << "[Atoms] (Angs)" << endl;
+  for (int i=0;i<natoms;i++)
+  {
+    string aname = get_aname(atno[i]);
+    double x1 = coords[3*i]*B2A; double y1 = coords[3*i+1]*B2A; double z1 = coords[3*i+2]*B2A;
+    outfile << " " << aname << "     " << i+1 << "   " << atno[i] << "   ";
+    outfile << x1 << "   " << y1 << "   " << z1 << endl;
+  }
+
+
+  outfile << "[STO]" << endl;
+
+  for (int i=0;i<natoms;i++)
+  {
+    double x1 = coords[3*i]; double y1 = coords[3*i+1]; double z1 = coords[3*i+2];
+
+    for (int j=0;j<N;j++)
+    {
+      if (close_val(x1,basis[j][5]) && close_val(y1,basis[j][6]) && close_val(z1,basis[j][7]))
+      {
+        int n1 = basis[j][0]; int l1 = basis[j][1]; int m1 = basis[j][2]; double zeta = basis[j][3]*B2A; double norm = basis[j][4];
+        int nx,ny,nz,nr;
+        get_nxyzr(n1,l1,m1,nx,ny,nz,nr);
+
+        if (l1<2)
+        {
+          outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+          outfile << nr << " " << zeta << " " << norm << " " << endl;
+        }
+        else if (l1==2)
+        {
+         //for 2 of the 5 3d functions, divide into parts
+          if (nz==2)
+          {
+           //3dz2 --> 2.z2 - x2 - y2
+            nz = 2; nx = ny = 0;
+            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+            outfile << nr << " " << zeta << " " << 2.*norm << " " << endl;
+            nz = 0; nx = 2; ny = 0;
+            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+            outfile << nr << " " << zeta << " " << norm << " " << endl;
+            nz = 0; nx = 0; ny = 2;
+            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+            outfile << nr << " " << zeta << " " << norm << " " << endl;
+          }
+          else if (nx==2 && ny==2)
+          {
+            //x2-y2 --> x2 - y2
+            nx = 2; ny = 0;
+            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+            outfile << nr << " " << zeta << " " << norm << " " << endl;
+            nx = 0; ny = 2;
+            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+            outfile << nr << " " << zeta << " " << norm << " " << endl;
+          }
+          else
+          {
+            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+            outfile << nr << " " << zeta << " " << norm << " " << endl;
+          }
+        } //if l1==2
+
+       #if 0
+       //just drop most f or higher functions
+        else if (l1==3)
+        {
+         //m==-2
+          if (nx==1 && ny==1 && nz==1) //fxyz is the simplest 4f orbital
+          {
+            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+            outfile << nr << " " << zeta << " " << norm << " " << endl;
+          }
+         //m==0
+          if (nz==3)
+          {
+           //5z3
+            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+            outfile << nr << " " << zeta << " " << 5.*norm << " " << endl;
+           //-3zr2
+            nz = 1; nr = 2;
+            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
+            outfile << nr << " " << zeta << " " << -3.*norm << " " << endl;
+          }
+          else
+          {
+            //incomplete
+          }
+        } //if l1==3
+       #endif
+
+      }
+    }
+  }
+
+  outfile << "[MO]" << endl;
+  for (int i=0;i<N;i++)
+  {
+    int occ = 0; if (i<No) occ = 1;
+
+    outfile << "Sym=X" << endl;
+    outfile << "Ene=" << eig[i] << endl;
+    outfile << "Spin=Alpha" << endl;
+    outfile << "Occup=" << occ << endl;
+
+    int wj = 0;
+    for (int j=0;j<N;j++)
+    {
+      int n1 = basis[j][0]; int l1 = basis[j][1]; int m1 = basis[j][2];
+      int nx,ny,nz,nr;
+      get_nxyzr(n1,l1,m1,nx,ny,nz,nr);
+
+      if (l1<2)
+        outfile << " " << 1+wj++ << "  " << jCA[j*N+i] << endl;
+      else if (l1==2)
+      {
+        if (nz==2)
+        {
+          outfile << " " << 1+wj++ << "  " <<  jCA[j*N+i] << endl;
+          outfile << " " << 1+wj++ << "  " << -jCA[j*N+i] << endl;
+          outfile << " " << 1+wj++ << "  " << -jCA[j*N+i] << endl;
+        }
+        else if (nx==2 && ny==2)
+        {
+          outfile << " " << 1+wj++ << "  " <<  jCA[j*N+i] << endl;
+          outfile << " " << 1+wj++ << "  " << -jCA[j*N+i] << endl;
+        }
+        else
+          outfile << " " << 1+wj++ << "  " << jCA[j*N+i] << endl;
+      }
+     #if 0
+      else if (l1==3)
+      {
+       //incomplete
+        if (nx==1 && ny==1 && nz==1) //m==-2
+          outfile << " " << 1+wj++ << "  " << jCA[j*N+i] << endl;
+        else if (nz==3) //m==0
+          outfile << " " << 1+wj++ << "  " << jCA[j*N+i] << endl;
+      }
+     #endif
+    }
+  }
+
+
+  outfile.close();
+
+  return;
+}
+
+void write_molden(bool gbasis, int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, string fname)
+{
+  int N = basis.size();
+  double eig[N];
+  for (int i=0;i<N;i++)
+    eig[i] = 0.;
+  return write_molden(gbasis,natoms,atno,coords,basis,jCA,No,eig,fname);
+}
+
+#if 0
+void write_molden_vcf(int natoms, int* atno, double* coords, vector<vector<vcf> > &vcfs, string fname)
+{
+  for (int n=0;n<natoms;n++)
+  for (int j=0;j<vcfs[n].size();j++)
+  if (vcfs[n][j].numerical)
+    return;
+
+  int nvcft = 0;
+  for (int n=0;n<natoms;n++)
+    nvcft += vcfs[n].size();
+  bool found_mu = 0;
+  for (int n=0;n<natoms;n++)
+  for (int j=0;j<vcfs[n].size();j++)
+  if (vcfs[n][j].mu)
+  { found_mu = 1; nvcft++; }
+  if (found_mu) printf(" WARNING: mu lobes from p orbitals in molden output \n");
+
+  int N = nvcft;
+  int N2 = N*N;
+
+  double* jCA = new double[N2]();
+  double R12 = 1.;
+  if (natoms==2) { double a12 = coords[0]-coords[3]; double b12 = coords[1]-coords[4]; double c12 = coords[2]-coords[5]; R12 = sqrt(a12*a12+b12*b12+c12*c12); }
+
+  vector<vector<double> > vbasis;
+  int s1 = 0;
+  for (int n=0;n<natoms;n++)
+  {
+    float A1 = coords[3*n+0]; float B1 = coords[3*n+1]; float C1 = coords[3*n+2];
+
+    vector<vcf> vcf1 = vcfs[n];
+    int nvcf = vcf1.size();
+
+    int jc = 0;
+    for (int j=0;j<nvcf;j++)
+    {
+      float v1 = vcf1[j].v1;
+
+      int n1 = vcf1[j].n; int l1 = vcf1[j].l; int m1 = vcf1[j].m;
+      vector<double> b1;
+      for (int m=0;m<10;m++) b1.push_back(0);
+      b1[5] = A1; b1[6] = B1; b1[7] = C1;
+      b1[8] = atno[n]; b1[9] = atno[n];
+
+      b1[0] = n1; b1[1] = l1; b1[2] = m1;
+      b1[3] = vcf1[j].zeta;
+      b1[4] = vcf1[j].norm;
+
+      vbasis.push_back(b1);
+
+      if (vcf1[j].mu>0)
+      {
+       //ftns on 2 atoms
+        vbasis.back()[3] = b1[3] = vcf1[j].zeta/R12;
+        if (n==0) { b1[5] = coords[3]; b1[6] = coords[4]; b1[7] = coords[5]; }
+        if (n==1) { b1[5] = coords[0]; b1[6] = coords[1]; b1[7] = coords[2]; }
+        vbasis.push_back(b1);
+
+        jCA[(s1+jc)*N+0] = v1;
+        jCA[(s1+jc+1)*N+0] = -v1;
+
+        jc++;
+      }
+      else if (vcf1[j].Vatom==0) //can't write Vatom
+      {
+        jCA[(s1+jc)*N+0] = v1;
+       //to visualize individual components (all but first component)
+        jCA[(s1+jc)*N+s1+jc] = v1;
+      }
+      jc++;
+
+    } //loop j for atom n
+    s1 += jc;
+  }
+
+  write_molden(0,natoms,atno,coords,vbasis,jCA,1,fname);
+
+  delete [] jCA;
+
+  return;
+}
+#endif
+
+#if 0
+void write_molden_vcf(int natoms, int* atno, float* coordsf, vector<vector<vcf> > &vcfs, string fname)
+{
+  double coords[3*natoms];
+  for (int j=0;j<3*natoms;j++)
+    coords[j] = coordsf[j];
+  return write_molden_vcf(natoms,atno,coords,vcfs,fname);
+}
+#endif
diff --git a/src/libio/read.cpp b/src/libio/read.cpp
index a3e1eca..00a90f7 100644
--- a/src/libio/read.cpp
+++ b/src/libio/read.cpp
@@ -64,9 +64,23 @@ vector<string> split1(const string &s, char delim)
   return tokens;
 }
 
-string read_basis_text(string aname)
+vector<string> split2(const string &s, char delim)
+{
+  stringstream ss(s);
+  string item;
+  vector<string> tokens1;
+  while ((bool)getline(ss, item, delim))
+    tokens1.push_back(item);
+
+  vector<string> tokens2;
+  for (const string& t : tokens1)
+    if (!t.empty()) tokens2.push_back(t);
+
+  return tokens2;
+}
+
+string read_basis_text(string filename, string aname)
 {
-  string filename = "basis";
   ifstream infile;
   infile.open(filename.c_str());
   if (!infile)
@@ -105,6 +119,12 @@ string read_basis_text(string aname)
   return text;
 }
 
+string read_basis_text(string aname)
+{
+  string filename = "basis";
+  return read_basis_text(filename,aname);
+}
+
 double read_float(string filename)
 {
   ifstream infile;
@@ -1853,18 +1873,19 @@ bool read_MOI_from_file(int Mm, int M3, double** MOI, string filename, int prl)
       wi++;
     }
   }
-  
+
   // Verify that the expected number of rows were read
-  if (wi == Mm) { 
-    if (prl > 1) printf("   found all lines of MOI \n"); 
+  if (wi==Mm)
+  {
+    if (prl > 1) printf("   found all lines of MOI \n");
   }
   else printf(" MOI missing lines \n");
-  
+
   infile.close();
   return true;
 }
 
-double nuclear_repulsion(int natoms, int* atno, double* coords)
+double nuclear_repulsion(int natoms, double* atno, double* coords)
 {
   double Enn = 0;
 
@@ -1881,6 +1902,19 @@ double nuclear_repulsion(int natoms, int* atno, double* coords)
   return Enn;
 }
 
+double nuclear_repulsion(int natoms, int* atno, double* coords)
+{
+  double* atnod = new double[natoms];
+  for (int n=0;n<natoms;n++)
+    atnod[n] = atno[n];
+
+  double Enn = nuclear_repulsion(natoms,atnod,coords);
+
+  delete [] atnod;
+
+  return Enn;
+}
+
 double nuclear_repulsion(int natoms, int* atno, float* coordsf)
 {
   double Enn = 0;
@@ -2547,10 +2581,10 @@ int read_input(string filename, bool gbasis, vector<vector<double> >& basis, vec
   while (!infile.eof())
   {
     (bool)getline(infile, line);
-    vector<string> tok_line_raw = split1(line, ' ');
-    vector<string> tok_line;
-    for (const string& t : tok_line_raw)
-      if (!t.empty()) tok_line.push_back(t);
+    vector<string> tok_line = split2(line, ' ');
+  //  vector<string> tok_line;
+  //  for (const string& t : tok_line_raw)
+  //    if (!t.empty()) tok_line.push_back(t);
 
     //cout << " 2READ: " << line << endl;
     if (tok_line.size()>3)
@@ -2655,6 +2689,92 @@ int read_input(string filename, bool gbasis, vector<vector<double> >& basis, vec
   return natoms;
 }
 
+int read_geoms(int& ngeom1, int natoms, vector<double*>& coords, int prl)
+{
+  string filename = "GEOMS";
+  int N3 = 3*natoms;
+
+  ifstream infile;
+  infile.open(filename.c_str());
+  if (!infile)
+  {
+    printf("  couldn't open file: %s \n",filename.c_str());
+    return 0;
+  }
+
+  string line;
+  bool found = (bool)getline(infile, line);
+  vector<string> tok_line = split2(line,' ');
+  if (tok_line.size()<1)
+  {
+    printf("  missing info in %s first line \n",filename.c_str());
+    return 0;
+  }
+
+  int ngeom = stoi(tok_line[0]);
+  if (tok_line.size()>1)
+    ngeom1 = stoi(tok_line[1]);
+  else
+    ngeom1 = ngeom;
+
+  if (prl>1)
+    printf(" looking for %i (%i) geometries with %i atoms \n",ngeom,ngeom1,natoms);
+  if (ngeom<1)
+  {
+    infile.close();
+    return 0;
+  }
+
+  {
+    for (int n=0;n<coords.size();n++)
+      delete [] coords[n];
+    coords.clear();
+  }
+
+  int wa = 0;
+  int wg = 0;
+  double* coords1 = new double[N3]();
+  while (!infile.eof())
+  {
+    (bool)getline(infile, line);
+    vector<string> tok_line = split2(line,' ');
+
+    if (tok_line.size()<4) { printf("\n ERROR: bad geometry \n"); exit(-1); }
+
+    coords1[3*wa+0] = stod(tok_line[1])*A2B;
+    coords1[3*wa+1] = stod(tok_line[2])*A2B;
+    coords1[3*wa+2] = stod(tok_line[3])*A2B;
+
+    wa++;
+    if (wa==natoms)
+    {
+      coords.push_back(coords1);
+      wa = 0; wg++;
+      if (wg==ngeom) break;
+
+      coords1 = new double[N3]();
+    }
+  }
+  if (wg<ngeom)
+    delete [] coords1;
+
+  infile.close();
+
+  if (prl>0)
+  {
+    printf(" geometries: \n");
+    for (int n=0;n<ngeom;n++)
+    {
+      for (int i=0;i<natoms;i++)
+        printf("  %8.5f %8.5f %8.5f \n",coords[n][3*i+0],coords[n][3*i+1],coords[n][3*i+2]);
+      printf("\n");
+    }
+  }
+
+  return ngeom;
+}
+
+
 void print_basis(int natoms, vector<vector<double> >& basis, vector<vector<double> >& basis_aux, int prl)
 {
   if (prl<0) return;
diff --git a/src/libio/write.cpp b/src/libio/write.cpp
index ff31bcc..093f741 100644
--- a/src/libio/write.cpp
+++ b/src/libio/write.cpp
@@ -9,6 +9,12 @@
 #define A2B 1.8897261
 
 
+bool close_val(double v1, double v2)
+{
+  if (fabs(v1-v2)<1.e-4) return true;
+  return false;
+}
+
 string get_iarray_name(short type1, short type2, short i1)
 {
   string ft;
@@ -72,7 +78,7 @@ void write_iarray(short type1, short type2, short i1, int s1, int s2, float* A)
 
 void write_iarray(short type1, short type2, short i1, int s1, int s2, double* A)
 {
-  printf("\n shouldn't be here (for now) \n"); 
+  printf("\n shouldn't be here (for now) \n");
 
   string filename = get_iarray_name(type1,type2,i1);
   //printf("  writing %8s \n",filename.c_str());
@@ -118,625 +124,6 @@ void save_geoms(int natoms, int* atno, vector<float> E, vector<float*> geoms, st
   return;
 }
 
-void get_nxyzr(int n1, int l1, int m1, int& nx, int& ny, int& nz, int& nr)
-{
-  nx=ny=nz=nr=0;
-  if (n1==1)
-  {
-  }
-  else if (n1==2)
-  {
-    if (l1==0)
-      nr = 1;
-    else if (l1==1)
-    {
-      if (m1==0)      nz = 1;
-      else if (m1==1) nx = 1;
-      else            ny = 1;
-    }
-  }
-  else if (n1==3)
-  {
-    if (l1==0)
-      nr = 2;
-    else if (l1==1)
-    {
-      nr = 1;
-      if (m1==0)      nz = 1;
-      else if (m1==1) nx = 1;
-      else            ny = 1;
-    }
-    else if (l1==2)
-    {
-      nr = 0;
-     #if CART_D
-      if (m1==0) { nx = 2; }
-      if (m1==1) { ny = 2; }
-      if (m1==2) { nz = 2; }
-      if (m1==3) { nx = 1; ny = 1; }
-      if (m1==4) { nx = 1; nz = 1; }
-      if (m1==5) { ny = 1; nz = 1; }
-     #else
-      if (m1==-2) { nx = 1; ny = 1; }
-      if (m1==-1) { ny = 1; nz = 1; }
-      if (m1== 0) { nz = 2; }
-      if (m1== 1) { nx = 1; nz = 1; }
-      if (m1== 2) { nx = 2; ny = 2; }
-     #endif
-    }
-  }
-  else if (n1==4)
-  {
-    //printf(" WARNING: n=4 MO print incomplete \n");
-    if (l1==0)
-      nr = 3;
-    else if (l1==1)
-    {
-      nr = 2;
-      if (m1==0)      nz = 1;
-      else if (m1==1) nx = 1;
-      else            ny = 1;
-    }
-    else if (l1==2)
-    {
-      nr = 1;
-     #if CART_D
-      if (m1==0) { nx = 2; }
-      if (m1==1) { ny = 2; }
-      if (m1==2) { nz = 2; }
-      if (m1==3) { nx = 1; ny = 1; }
-      if (m1==4) { nx = 1; nz = 1; }
-      if (m1==5) { ny = 1; nz = 1; }
-     #else
-      if (m1==-2) { nx = 1; ny = 1; }
-      if (m1==-1) { ny = 1; nz = 1; }
-      if (m1== 0) { nz = 2; }
-      if (m1== 1) { nx = 1; nz = 1; }
-      if (m1== 2) { nx = 2; ny = 2; }
-     #endif
-    }
-    else if (l1==3)
-    {
-     //INCOMPLETE
-      if (m1==-2) { nx = 1; ny = 1; nz = 1; }
-      else if (m1==0) { nz = 3; }
-      else
-        nx = nx = nz = 444;
-    }
-
-  }
-  return;
-}
-
-void write_molden_g(int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, double* eig, string fname)
-{
-  int N = basis.size();
-
-  string filename = fname;
-  ofstream outfile;
-  outfile.open(filename.c_str());
-
-  double B2A = 1./A2B;
-
-  outfile << fixed << setprecision(8);
-
-  outfile << "[Molden Format]" << endl;
-  outfile << "[Atoms] (Angs)" << endl;
-  for (int i=0;i<natoms;i++)
-  {
-    string aname = get_aname(atno[i]);
-    double x1 = coords[3*i]*B2A; double y1 = coords[3*i+1]*B2A; double z1 = coords[3*i+2]*B2A;
-    outfile << " " << aname << "     " << i+1 << "   " << atno[i] << "   ";
-    outfile << x1 << "   " << y1 << "   " << z1 << endl;
-  }
-
-
-  outfile << "[GTO]" << endl;
-
-  for (int i=0;i<natoms;i++)
-  {
-    //double x1 = coords[3*i]; double y1 = coords[3*i+1]; double z1 = coords[3*i+2];
-
-    string aname = get_aname(atno[i]);
-    string b1 = read_basis_text(aname);
-    outfile << i+1 << "    0 " << endl;
-    outfile << b1 << endl;
-      //outfile << nr << " " << zeta << " " << norm << " " << endl;
-  }
-
-  outfile << "[MO]" << endl;
-  for (int i=0;i<N;i++)
-  {
-    int occ = 0; if (i<No) occ = 1;
-
-    outfile << "Sym=X" << endl;
-    outfile << "Ene=" << eig[i] << endl;
-    outfile << "Spin=Alpha" << endl;
-    outfile << "Occup=" << occ << endl;
-
-    for (int j=0;j<N;j++)
-      outfile << " " << j+1 << "  " << jCA[j*N+i] << endl;
-  }
-
-  outfile << " [5D7F] " << endl;
-
-
-  outfile.close();
-
-  return;
-}
-
-bool close_val(double v1, double v2)
-{
-  if (fabs(v1-v2)<1.e-4) return true;
-  return false;
-}
-
-void write_molden_ss(int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, string fname)
-{
-  int N = basis.size();
-
-  string filename = fname;
-  ofstream outfile;
-  outfile.open(filename.c_str());
-
-  double B2A = 1./A2B;
-
-  outfile << fixed << setprecision(8);
-
-  outfile << "[Molden Format]" << endl;
-  outfile << "[Atoms] (Angs)" << endl;
-  for (int i=0;i<natoms;i++)
-  {
-    string aname = get_aname(atno[i]);
-    double x1 = coords[3*i]*B2A; double y1 = coords[3*i+1]*B2A; double z1 = coords[3*i+2]*B2A;
-    outfile << " " << aname << "     " << i+1 << "   " << atno[i] << "   ";
-    outfile << x1 << "   " << y1 << "   " << z1 << endl;
-  }
-
-
-  outfile << "[STO]" << endl;
-
-  int nss = 4;
-  for (int i=0;i<natoms;i++)
-  {
-    double x1 = coords[3*i]; double y1 = coords[3*i+1]; double z1 = coords[3*i+2];
-
-    for (int j=0;j<N;j++)
-    {
-      if (close_val(x1,basis[j][5]) && close_val(y1,basis[j][6]) && close_val(z1,basis[j][7]))
-      {
-        int n1 = basis[j][0]; int l1 = basis[j][1]; int m1 = basis[j][2]; double zeta = basis[j][3]*B2A; double norm = basis[j][4];
-        int nx,ny,nz,nr;
-        get_nxyzr(n1,l1,m1,nx,ny,nz,nr);
-
-        if (l1<2)
-        {
-          for (int n=0;n<nss;n++)
-          {
-           //nr+n-l1?
-            double norm1 = norm*pow(zeta,n); if (n==2) norm1 *= 0.4; if (n==3) norm1 /= 15.;
-            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-            outfile << nr+n << " " << zeta << " " << norm1 << " " << endl;
-          }
-        }
-        else if (l1==2)
-        {
-         //for 2 of the 5 3d functions, divide into parts
-          if (nz==2)
-          {
-           //3dz2 --> 2.z2 - x2 - y2
-            for (int n=0;n<nss;n++)
-            {
-              double norm1 = norm*pow(zeta,n); if (n==2) norm1 *= 0.4; if (n==3) norm1 /= 15.;
-              nz = 2; nx = ny = 0;
-              outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-              outfile << nr+n << " " << zeta << " " << 2.*norm1 << " " << endl;
-              nz = 0; nx = 2; ny = 0;
-              outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-              outfile << nr+n << " " << zeta << " " << norm1 << " " << endl;
-              nz = 0; nx = 0; ny = 2;
-              outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-              outfile << nr+n << " " << zeta << " " << norm1 << " " << endl;
-            }
-          }
-          else if (nx==2 && ny==2)
-          {
-            //x2-y2 --> x2 - y2
-            for (int n=0;n<nss;n++)
-            {
-              double norm1 = norm*pow(zeta,n); if (n==2) norm1 *= 0.4; if (n==3) norm1 /= 15.;
-              nx = 2; ny = 0;
-              outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-              outfile << nr+n << " " << zeta << " " << norm1 << " " << endl;
-              nx = 0; ny = 2;
-              outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-              outfile << nr+n << " " << zeta << " " << norm1 << " " << endl;
-            }
-          }
-          else
-          {
-            for (int n=0;n<nss;n++)
-            {
-              double norm1 = norm*pow(zeta,n); if (n==2) norm1 *= 0.4; if (n==3) norm1 /= 15.;
-              outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-              outfile << nr+n << " " << zeta << " " << norm1 << " " << endl;
-            }
-          }
-        } //if l1==2
-
-       #if 0
-       //just drop most f or higher functions
-        else if (l1==3)
-        {
-         //m==-2
-          if (nx==1 && ny==1 && nz==1) //fxyz is the simplest 4f orbital
-          {
-            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-            outfile << nr << " " << zeta << " " << norm << " " << endl;
-          }
-         //m==0
-          if (nz==3)
-          {
-           //5z3
-            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-            outfile << nr << " " << zeta << " " << 5.*norm << " " << endl;
-           //-3zr2
-            nz = 1; nr = 2;
-            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-            outfile << nr << " " << zeta << " " << -3.*norm << " " << endl;
-          }
-          else
-          {
-            //incomplete
-          }
-        } //if l1==3
-       #endif
-
-      }
-    }
-  }
-
-  outfile << "[MO]" << endl;
-  for (int i=0;i<N;i++)
-  {
-    int occ = 0; if (i<No) occ = 1;
-
-    outfile << "Sym=X" << endl;
-    outfile << "Ene=0.0" << endl;
-    outfile << "Spin=Alpha" << endl;
-    outfile << "Occup=" << occ << endl;
-
-    int wj = 0;
-    for (int j=0;j<N;j++)
-    {
-      int n1 = basis[j][0]; int l1 = basis[j][1]; int m1 = basis[j][2];
-      int nx,ny,nz,nr;
-      get_nxyzr(n1,l1,m1,nx,ny,nz,nr);
-
-      if (l1<2)
-      {
-        for (int n=0;n<nss;n++)
-          outfile << " " << 1+wj++ << "  " << jCA[j*N+i] << endl;
-      }
-      else if (l1==2)
-      {
-        if (nz==2)
-        {
-          for (int n=0;n<nss;n++)
-          {
-            outfile << " " << 1+wj++ << "  " <<  jCA[j*N+i] << endl;
-            outfile << " " << 1+wj++ << "  " << -jCA[j*N+i] << endl;
-            outfile << " " << 1+wj++ << "  " << -jCA[j*N+i] << endl;
-          }
-        }
-        else if (nx==2 && ny==2)
-        {
-          for (int n=0;n<nss;n++)
-          {
-            outfile << " " << 1+wj++ << "  " <<  jCA[j*N+i] << endl;
-            outfile << " " << 1+wj++ << "  " << -jCA[j*N+i] << endl;
-          }
-        }
-        else
-        {
-          for (int n=0;n<nss;n++)
-            outfile << " " << 1+wj++ << "  " << jCA[j*N+i] << endl;
-        }
-      }
-     #if 0
-      else if (l1==3)
-      {
-       //incomplete
-        if (nx==1 && ny==1 && nz==1) //m==-2
-          outfile << " " << 1+wj++ << "  " << jCA[j*N+i] << endl;
-        else if (nz==3) //m==0
-          outfile << " " << 1+wj++ << "  " << jCA[j*N+i] << endl;
-      }
-     #endif
-    }
-  }
-
-  outfile.close();
-
-  return;
-}
-
-void write_molden(bool gbasis, int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, double* eig, string fname)
-{
-  if (gbasis) return write_molden_g(natoms,atno,coords,basis,jCA,No,eig,fname);
-  if (basis[0].size()>10) return write_molden_ss(natoms,atno,coords,basis,jCA,No,fname);
-
-  int N = basis.size();
-
-  bool missing_ftns = 0;
-  for (int j=0;j<N;j++)
-  if (basis[j][0]>3)
-    missing_ftns = 1;
-  if (missing_ftns)
-    printf(" WARNING: n=4 MO print incomplete \n");
-  //printf("\n WARNING: no f function printing to molden \n");
-
-  string filename = fname;
-  ofstream outfile;
-  outfile.open(filename.c_str());
-
-  double B2A = 1./A2B;
-
-  outfile << fixed << setprecision(8);
-
-  outfile << "[Molden Format]" << endl;
-  outfile << "[Atoms] (Angs)" << endl;
-  for (int i=0;i<natoms;i++)
-  {
-    string aname = get_aname(atno[i]);
-    double x1 = coords[3*i]*B2A; double y1 = coords[3*i+1]*B2A; double z1 = coords[3*i+2]*B2A;
-    outfile << " " << aname << "     " << i+1 << "   " << atno[i] << "   ";
-    outfile << x1 << "   " << y1 << "   " << z1 << endl;
-  }
-
-
-  outfile << "[STO]" << endl;
-
-  for (int i=0;i<natoms;i++)
-  {
-    double x1 = coords[3*i]; double y1 = coords[3*i+1]; double z1 = coords[3*i+2];
-
-    for (int j=0;j<N;j++)
-    {
-      if (close_val(x1,basis[j][5]) && close_val(y1,basis[j][6]) && close_val(z1,basis[j][7]))
-      {
-        int n1 = basis[j][0]; int l1 = basis[j][1]; int m1 = basis[j][2]; double zeta = basis[j][3]*B2A; double norm = basis[j][4];
-        int nx,ny,nz,nr;
-        get_nxyzr(n1,l1,m1,nx,ny,nz,nr);
-
-        if (l1<2)
-        {
-          outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-          outfile << nr << " " << zeta << " " << norm << " " << endl;
-        }
-        else if (l1==2)
-        {
-         //for 2 of the 5 3d functions, divide into parts
-          if (nz==2)
-          {
-           //3dz2 --> 2.z2 - x2 - y2
-            nz = 2; nx = ny = 0;
-            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-            outfile << nr << " " << zeta << " " << 2.*norm << " " << endl;
-            nz = 0; nx = 2; ny = 0;
-            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-            outfile << nr << " " << zeta << " " << norm << " " << endl;
-            nz = 0; nx = 0; ny = 2;
-            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-            outfile << nr << " " << zeta << " " << norm << " " << endl;
-          }
-          else if (nx==2 && ny==2)
-          {
-            //x2-y2 --> x2 - y2
-            nx = 2; ny = 0;
-            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-            outfile << nr << " " << zeta << " " << norm << " " << endl;
-            nx = 0; ny = 2;
-            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-            outfile << nr << " " << zeta << " " << norm << " " << endl;
-          }
-          else
-          {
-            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-            outfile << nr << " " << zeta << " " << norm << " " << endl;
-          }
-        } //if l1==2
-
-       #if 0
-       //just drop most f or higher functions
-        else if (l1==3)
-        {
-         //m==-2
-          if (nx==1 && ny==1 && nz==1) //fxyz is the simplest 4f orbital
-          {
-            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-            outfile << nr << " " << zeta << " " << norm << " " << endl;
-          }
-         //m==0
-          if (nz==3)
-          {
-           //5z3
-            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-            outfile << nr << " " << zeta << " " << 5.*norm << " " << endl;
-           //-3zr2
-            nz = 1; nr = 2;
-            outfile << i+1 << " " << nx << " " << ny << " " << nz << " ";
-            outfile << nr << " " << zeta << " " << -3.*norm << " " << endl;
-          }
-          else
-          {
-            //incomplete
-          }
-        } //if l1==3
-       #endif
-
-      }
-    }
-  }
-
-  outfile << "[MO]" << endl;
-  for (int i=0;i<N;i++)
-  {
-    int occ = 0; if (i<No) occ = 1;
-
-    outfile << "Sym=X" << endl;
-    outfile << "Ene=" << eig[i] << endl;
-    outfile << "Spin=Alpha" << endl;
-    outfile << "Occup=" << occ << endl;
-
-    int wj = 0;
-    for (int j=0;j<N;j++)
-    {
-      int n1 = basis[j][0]; int l1 = basis[j][1]; int m1 = basis[j][2];
-      int nx,ny,nz,nr;
-      get_nxyzr(n1,l1,m1,nx,ny,nz,nr);
-
-      if (l1<2)
-        outfile << " " << 1+wj++ << "  " << jCA[j*N+i] << endl;
-      else if (l1==2)
-      {
-        if (nz==2)
-        {
-          outfile << " " << 1+wj++ << "  " <<  jCA[j*N+i] << endl;
-          outfile << " " << 1+wj++ << "  " << -jCA[j*N+i] << endl;
-          outfile << " " << 1+wj++ << "  " << -jCA[j*N+i] << endl;
-        }
-        else if (nx==2 && ny==2)
-        {
-          outfile << " " << 1+wj++ << "  " <<  jCA[j*N+i] << endl;
-          outfile << " " << 1+wj++ << "  " << -jCA[j*N+i] << endl;
-        }
-        else
-          outfile << " " << 1+wj++ << "  " << jCA[j*N+i] << endl;
-      }
-     #if 0
-      else if (l1==3)
-      {
-       //incomplete
-        if (nx==1 && ny==1 && nz==1) //m==-2
-          outfile << " " << 1+wj++ << "  " << jCA[j*N+i] << endl;
-        else if (nz==3) //m==0
-          outfile << " " << 1+wj++ << "  " << jCA[j*N+i] << endl;
-      }
-     #endif
-    }
-  }
-
-
-  outfile.close();
-
-  return;
-}
-
-void write_molden(bool gbasis, int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, string fname)
-{
-  int N = basis.size();
-  double eig[N];
-  for (int i=0;i<N;i++)
-    eig[i] = 0.;
-  return write_molden(gbasis,natoms,atno,coords,basis,jCA,No,eig,fname);
-}
-
-#if 0
-void write_molden_vcf(int natoms, int* atno, double* coords, vector<vector<vcf> > &vcfs, string fname)
-{
-  for (int n=0;n<natoms;n++)
-  for (int j=0;j<vcfs[n].size();j++)
-  if (vcfs[n][j].numerical)
-    return;
-
-  int nvcft = 0;
-  for (int n=0;n<natoms;n++)
-    nvcft += vcfs[n].size();
-  bool found_mu = 0;
-  for (int n=0;n<natoms;n++)
-  for (int j=0;j<vcfs[n].size();j++)
-  if (vcfs[n][j].mu)
-  { found_mu = 1; nvcft++; }
-  if (found_mu) printf(" WARNING: mu lobes from p orbitals in molden output \n");
-
-  int N = nvcft;
-  int N2 = N*N;
-
-  double* jCA = new double[N2]();
-  double R12 = 1.;
-  if (natoms==2) { double a12 = coords[0]-coords[3]; double b12 = coords[1]-coords[4]; double c12 = coords[2]-coords[5]; R12 = sqrt(a12*a12+b12*b12+c12*c12); }
-
-  vector<vector<double> > vbasis;
-  int s1 = 0;
-  for (int n=0;n<natoms;n++)
-  {
-    float A1 = coords[3*n+0]; float B1 = coords[3*n+1]; float C1 = coords[3*n+2];
-
-    vector<vcf> vcf1 = vcfs[n];
-    int nvcf = vcf1.size();
-
-    int jc = 0;
-    for (int j=0;j<nvcf;j++)
-    {
-      float v1 = vcf1[j].v1;
-
-      int n1 = vcf1[j].n; int l1 = vcf1[j].l; int m1 = vcf1[j].m;
-      vector<double> b1;
-      for (int m=0;m<10;m++) b1.push_back(0);
-      b1[5] = A1; b1[6] = B1; b1[7] = C1;
-      b1[8] = atno[n]; b1[9] = atno[n];
-
-      b1[0] = n1; b1[1] = l1; b1[2] = m1;
-      b1[3] = vcf1[j].zeta;
-      b1[4] = vcf1[j].norm;
-
-      vbasis.push_back(b1);
-
-      if (vcf1[j].mu>0)
-      {
-       //ftns on 2 atoms
-        vbasis.back()[3] = b1[3] = vcf1[j].zeta/R12;
-        if (n==0) { b1[5] = coords[3]; b1[6] = coords[4]; b1[7] = coords[5]; }
-        if (n==1) { b1[5] = coords[0]; b1[6] = coords[1]; b1[7] = coords[2]; }
-        vbasis.push_back(b1);
-
-        jCA[(s1+jc)*N+0] = v1;
-        jCA[(s1+jc+1)*N+0] = -v1;
-
-        jc++;
-      }
-      else if (vcf1[j].Vatom==0) //can't write Vatom
-      {
-        jCA[(s1+jc)*N+0] = v1;
-       //to visualize individual components (all but first component)
-        jCA[(s1+jc)*N+s1+jc] = v1;
-      }
-      jc++;
-
-    } //loop j for atom n
-    s1 += jc;
-  }
-
-  write_molden(0,natoms,atno,coords,vbasis,jCA,1,fname);
-
-  delete [] jCA;
-
-  return;
-}
-#endif
-
-#if 0
-void write_molden_vcf(int natoms, int* atno, float* coordsf, vector<vector<vcf> > &vcfs, string fname)
-{
-  double coords[3*natoms];
-  for (int j=0;j<3*natoms;j++)
-    coords[j] = coordsf[j];
-  return write_molden_vcf(natoms,atno,coords,vcfs,fname);
-}
-#endif
-
 void write_graph(int size, double* h, string filename)
 {
   ofstream outfile;
@@ -2110,4 +1497,4 @@ void write_mo_grid(int natoms, int* atno, double* coords, int nrad, int gsa, vec
   #pragma acc exit data delete(gridf[0:gs6],wtf[0:gsa])
   delete [] gridf;
   delete [] wtf;
-}
\ No newline at end of file
+}