diff --git a/ChangeLog b/ChangeLog
index c37675df..aae0a35e 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -3,6 +3,15 @@
 -Name
 -changes
 
+--------------
+Jan 27, 2023
+Name: Xin Jing
+Changes: (initialization.c)
+1. Fix a bug when the cell is a rotated cuboidal, use non-orthogonal calculation for that
+2. Add the warning for user to change the cell to get best performance
+3. Add warnings into .out file
+4. Remove SCF scriteria for QE and add warnings.
+
 --------------
 Oct 31, 2023
 Name: Boqin Zhang
diff --git a/src/electronicGroundState.c b/src/electronicGroundState.c
index 0d473f7b..3dd7f451 100644
--- a/src/electronicGroundState.c
+++ b/src/electronicGroundState.c
@@ -636,12 +636,6 @@ void scf_loop(SPARC_OBJ *pSPARC) {
         }
 #endif
         
-        // used for QE scf error, save input rho_in and phi_in
-        if (pSPARC->scf_err_type == 1) {
-            memcpy(pSPARC->phi_dmcomm_phi_in, pSPARC->elecstPotential, DMnd * sizeof(double));
-            memcpy(pSPARC->rho_dmcomm_phi_in, pSPARC->electronDens   , DMnd * sizeof(double));
-		}
-
         // update Veff_loc_dmcomm_phi_in
         if (pSPARC->MixingVariable == 1) {
             double *Veff_out = (pSPARC->spin_typ == 2) ? pSPARC->Veff_dia_loc_dmcomm_phi : pSPARC->Veff_loc_dmcomm_phi;
@@ -770,12 +764,8 @@ void scf_loop(SPARC_OBJ *pSPARC) {
         int scf_conv = 0;
 
         // find SCF error
-        if (pSPARC->scf_err_type == 0) { // default
-            Evaluate_scf_error(pSPARC, &error, &scf_conv);
-        } else if (pSPARC->scf_err_type == 1) { // QE scf err: conv_thr
-            Evaluate_QE_scf_error(pSPARC, &error, &scf_conv);
-        }
-        
+        Evaluate_scf_error(pSPARC, &error, &scf_conv);
+
         // check if Etot is NaN
         if (pSPARC->Etot != pSPARC->Etot) {
             if (!rank) printf("ERROR: Etot is NaN\n");
@@ -871,10 +861,6 @@ void scf_loop(SPARC_OBJ *pSPARC) {
             exit(EXIT_FAILURE);
         }
         fprintf(output_fp,"Total number of SCF: %-6d\n",SCFcount);
-        // for density mixing, extra poisson solve is needed
-        if (pSPARC->scf_err_type == 1 && pSPARC->MixingVariable == 0) { 
-            fprintf(output_fp,"Extra time for evaluating QE SCF Error: %.3f (sec)\n", pSPARC->t_qe_extra);
-        }
         fclose(output_fp);
     }
 
@@ -1030,55 +1016,6 @@ void Evaluate_scf_error(SPARC_OBJ *pSPARC, double *scf_error, int *scf_conv) {
 }
 
 
-
-/**
- * @brief Evaluate the scf error defined in Quantum Espresso.
- *
- *        Find the scf error defined in Quantum Espresso. QE implements 
- *        Eq.(A.7) of the reference paper, with a slight modification: 
- *          conv_thr = 4 \pi e^2 \Omega \sum_G |\Delta \rho(G)|^2 / G^2
- *        This is equivalent to 2 * Eq.(A.6) in the reference paper.
- *
- * @ref   P Giannozzi et al, J. Phys.:Condens. Matter 21(2009) 395502.
- */
-void Evaluate_QE_scf_error(SPARC_OBJ *pSPARC, double *scf_error, int *scf_conv) 
-{
-    int rank;
-    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
-
-    // update phi_out for density mixing
-    if (pSPARC->MixingVariable == 0) { // desity mixing
-        double t1, t2;
-        t1 = MPI_Wtime();
-        // solve the poisson equation for electrostatic potential, "phi"
-        Calculate_elecstPotential(pSPARC);
-        t2 = MPI_Wtime();
-        if (!rank) printf("QE scf error: update phi_out took %.3f ms\n", (t2-t1)*1e3); 
-        pSPARC->t_qe_extra += (t2 - t1);
-    }
-
-    double error = 0.0;
-    if (pSPARC->dmcomm_phi != MPI_COMM_NULL) {
-        int i;
-        int DMnd = pSPARC->Nd_d;
-        double loc_err = 0.0;
-        for (i = 0; i < DMnd; i++) {
-            loc_err += (pSPARC->electronDens[i]    - pSPARC->rho_dmcomm_phi_in[i]) * 
-                       (pSPARC->elecstPotential[i] - pSPARC->phi_dmcomm_phi_in[i]);
-        }
-        loc_err = fabs(loc_err * pSPARC->dV); // in case error is not numerically positive
-        MPI_Allreduce(&loc_err, &error, 1, MPI_DOUBLE, MPI_SUM, pSPARC->dmcomm_phi);
-    }
-
-    MPI_Bcast(&error, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);   
-    // output
-    *scf_error = error;
-    *scf_conv  = (pSPARC->usefock % 2 == 1) 
-                ? ((int) (error < pSPARC->TOL_SCF_INIT)) : ((int) (error < pSPARC->TOL_SCF));
-}
-
-
-
 /**
 * @ brief find net magnetization of the system
 **/
diff --git a/src/electrostatics.c b/src/electrostatics.c
index 45cd9547..ecd5d9c8 100644
--- a/src/electrostatics.c
+++ b/src/electrostatics.c
@@ -586,7 +586,8 @@ void GetInfluencingAtoms(SPARC_OBJ *pSPARC) {
     double Lx, Ly, Lz, rbx, rby, rbz, x0, y0, z0, x0_i, y0_i, z0_i;
     double DMxs, DMxe, DMys, DMye, DMzs, DMze;
     int rb_xl, rb_xr, rb_yl, rb_yr, rb_zl, rb_zr;
-    int isRbOut[3] = {0,0,0}; // check rb-region of atom is out of the domain
+    int *isRbOut = pSPARC->isRbOut; // check rb-region of atom is out of the domain
+    isRbOut[0] = isRbOut[1] = isRbOut[2] = 0; 
     MPI_Comm grid_comm = pSPARC->dmcomm_phi;
     MPI_Comm_size(grid_comm, &nproc);
     MPI_Comm_rank(grid_comm, &rank);
@@ -757,9 +758,9 @@ void GetInfluencingAtoms(SPARC_OBJ *pSPARC) {
     } 
     if (pSPARC->BCx == 1 || pSPARC->BCy == 1 || pSPARC->BCz == 1) {        
         MPI_Allreduce(MPI_IN_PLACE, isRbOut, 3, MPI_INT, MPI_LAND, pSPARC->dmcomm_phi);
-        #ifdef DEBUG
-        if (!rank) printf("Is Rb region out? isRbOut = [%d, %d, %d]\n", isRbOut[0], isRbOut[1], isRbOut[2]);
-        #endif
+        if (isRbOut[0] || isRbOut[1] || isRbOut[2]) {
+            if (!rank) printf("WARNING: Atoms are too close to boundary for pseudocharge calculation.\n");
+        }
     }
 
 }
diff --git a/src/finalization.c b/src/finalization.c
index 220ef795..8f6aa4a8 100644
--- a/src/finalization.c
+++ b/src/finalization.c
@@ -220,12 +220,6 @@ void Free_basic(SPARC_OBJ *pSPARC) {
         if (pSPARC->MixingVariable == 0 && pSPARC->spin_typ) {
             free(pSPARC->electronDens_in);
         }
-
-        // for using QE scf error definition
-        if (pSPARC->scf_err_type == 1) {
-            free(pSPARC->rho_dmcomm_phi_in);
-            free(pSPARC->phi_dmcomm_phi_in);
-        }
         
         free(pSPARC->mixing_hist_Pfk);
     }
diff --git a/src/highT/sqParallelization.c b/src/highT/sqParallelization.c
index 5a54c094..f1231216 100644
--- a/src/highT/sqParallelization.c
+++ b/src/highT/sqParallelization.c
@@ -337,14 +337,6 @@ void Setup_Comms_SQ(SPARC_OBJ *pSPARC) {
             assert(pSPARC->Veff_loc_dmcomm_phi_in != NULL);
         }
 
-        // The following rho_in and phi_in are only used for evaluating QE scf errors
-        if (pSPARC->scf_err_type == 1) {
-            pSPARC->rho_dmcomm_phi_in = (double *)malloc(DMnd * sizeof(double));
-            assert(pSPARC->rho_dmcomm_phi_in != NULL);
-            pSPARC->phi_dmcomm_phi_in = (double *)malloc(DMnd * sizeof(double));
-            assert(pSPARC->phi_dmcomm_phi_in != NULL);
-        }
-
         // initialize electrostatic potential as random guess vector
         if (pSPARC->FixRandSeed == 1) {
             SeededRandVec(pSPARC->elecstPotential, pSPARC->DMVertices, gridsizes, -1.0, 1.0, 0);
diff --git a/src/include/initialization.h b/src/include/initialization.h
index 03e6f2e8..c542c78b 100644
--- a/src/include/initialization.h
+++ b/src/include/initialization.h
@@ -94,11 +94,6 @@ void Calculate_SplineDerivRadFun(SPARC_OBJ *pSPARC);
 
 void Cart2nonCart_transformMat(SPARC_OBJ *pSPARC);
 
-/**
- * @brief   Write the initialized parameters into the output file.
- *
- * @param pSPARC    The pointer that points to SPARC_OBJ type structure.
- */
  
 /*
  * @brief   function to convert non cartesian to cartesian coordinates
@@ -135,6 +130,8 @@ void CalculateDistance(SPARC_OBJ *pSPARC, double x, double y, double z, double x
 void write_output_init(SPARC_OBJ *pSPARC);
 
 
+void print_orthogonal_warning(SPARC_OBJ *pSPARC, FILE *output_fp);
+
 /**
  * @brief   Calculate k-points and the associated weights.
  *
diff --git a/src/include/isddft.h b/src/include/isddft.h
index b887cb0c..88966eec 100644
--- a/src/include/isddft.h
+++ b/src/include/isddft.h
@@ -222,6 +222,7 @@ typedef struct _SPARC_OBJ{
     int npNdx_kptcomm;  // number of processes in x-dir for creating Cartesian topology in kptcomm 
     int npNdy_kptcomm;  // number of processes in y-dir for creating Cartesian topology in kptcomm 
     int npNdz_kptcomm;  // number of processes in z-dir for creating Cartesian topology in kptcomm
+    int useDefaultParalFlag; // Flag for using default parallelization
     int FixRandSeed;    // flag to fix the random number seeds so that all random numbers generated in parallel 
                         // under MPI are the same as those generated in sequential execution
                         
@@ -346,6 +347,7 @@ typedef struct _SPARC_OBJ{
     int Nd_d_dmcomm;          // total number of grids of distributed domain in each dmcomm process (LOCAL)
     
     ATOM_INFLUENCE_OBJ *Atom_Influence_local; // atom info. for atoms that have local influence on the distributed domain (LOCAL)
+    int isRbOut[3];           // flag for if Rb region is outside domain for Dirichlet BC
     
     /* nonlocal */
     ATOM_NLOC_INFLUENCE_OBJ *Atom_Influence_nloc; // atom info. for atoms that have nonlocal influence on the distributed domain (LOCAL)
@@ -447,12 +449,7 @@ typedef struct _SPARC_OBJ{
     double *mixing_hist_Xk;      // residual matrix of Veff_loc, for mixing (LOCAL)
     double *mixing_hist_Fk;      // residual matrix of the residual of Veff_loc (LOCAL)
     double *mixing_hist_Pfk;     // the preconditioned residual distributed in phi-domain (LOCAL)
-
-    int    scf_err_type;        // scf error definition type
-    double t_qe_extra;          // // this is the extra unnecessary time we spent in order to evaluate QE scf error
-    // these two arrays are used only for evaluating QE scf error
-    double *rho_dmcomm_phi_in;  // input electron density distributed in phi-domain (LOCAL)
-    double *phi_dmcomm_phi_in;  // input electrostatic potential distributed in phi-domain (LOCAL)
+    
     double *psdChrgDens;          // pseudocharge density, "b" (LOCAL)
     double *psdChrgDens_ref;      // reference pseudocharge density, "b_ref" (LOCAL)
     double *Vc;                   // difference between reference pseudopotential V_ref and pseudopotential V, Vc = V_ref - V (LOCAL)
@@ -1058,7 +1055,6 @@ typedef struct _SPARC_INPUT_OBJ{
                         // under MPI are the same as those generated in sequential execution
                         // 0 - off (default), 1 - on
     int accuracy_level; // accuracy level, 1 - 'low', 2 - 'medium', 3 - 'high', 4 - 'extreme'
-    int scf_err_type;   // scf error definition type
     int MAXIT_SCF;      // max number of SCF iterations
     int MINIT_SCF;      // min number of SCF iterations
     int MAXIT_POISSON;  // max number of iterations for Poisson solver
diff --git a/src/initialization.c b/src/initialization.c
index f07a9140..8462b340 100644
--- a/src/initialization.c
+++ b/src/initialization.c
@@ -44,7 +44,7 @@
 #define min(x,y) ((x)<(y)?(x):(y))
 #define max(x,y) ((x)>(y)?(x):(y))
 
-#define N_MEMBR 175
+#define N_MEMBR 174
 
 
 /**
@@ -584,8 +584,7 @@ void set_defaults(SPARC_INPUT_OBJ *pSPARC_Input, SPARC_OBJ *pSPARC) {
     pSPARC_Input->Poisson_solver = 0;          // default AAR solver
     /* Iterations: tolerances and max_iters */
     pSPARC_Input->FixRandSeed = 0;            // default flag for fixing random numbers for MPI paralllelization
-    pSPARC_Input->accuracy_level = -1;        // default accuracy level (2 - 'medium', 1e-3 in energy and forces)
-    pSPARC_Input->scf_err_type = 0;           // default scf error definition: relative error in rho or veff
+    pSPARC_Input->accuracy_level = -1;        // default accuracy level (2 - 'medium', 1e-3 in energy and forces)    
     pSPARC_Input->MAXIT_SCF = 100;            // default maximum number of SCF iterations
     pSPARC_Input->MINIT_SCF = 2;              // default minimum number of SCF iterations
     pSPARC_Input->MAXIT_POISSON = 3000;       // default maximum number of iterations for Poisson solver
@@ -1179,7 +1178,6 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {
     pSPARC->Nchefsi = pSPARC_Input->Nchefsi;
     pSPARC->FixRandSeed = pSPARC_Input->FixRandSeed;
     pSPARC->accuracy_level = pSPARC_Input->accuracy_level;
-    pSPARC->scf_err_type = pSPARC_Input->scf_err_type;
     pSPARC->MAXIT_SCF = pSPARC_Input->MAXIT_SCF;
     pSPARC->MINIT_SCF = pSPARC_Input->MINIT_SCF;
     pSPARC->MAXIT_POISSON = pSPARC_Input->MAXIT_POISSON;
@@ -1502,7 +1500,7 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {
         printf("\nChecking existence of (%d) out file(s) took %.3f ms\n", i, (t2-t1)*1000);
 #endif
         if (i >= (int)(MAX_OUTPUT*0.75) && i <= MAX_OUTPUT) {
-            printf("\n#WARNING: There's a limit on total number of output files allowed!\n"
+            printf("\nWARNING: There's a limit on total number of output files allowed!\n"
                      "         After the limit is reached, SPARC will start using  the\n"
                      "         name provided exactly (without attached index) and old\n"
                      "         results will be overwritten! Please move some existing\n"
@@ -1512,7 +1510,7 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {
         }
 
         if (i > MAX_OUTPUT) {
-            printf("\n#WARNING: The limit of total number of output files is reached! \n"
+            printf("\nWARNING: The limit of total number of output files is reached! \n"
                      "         Current output name (without suffix): %s\n\n", pSPARC->filename_out);
         } else if (i > 0) {
             char tempchar[L_STRING];
@@ -1609,6 +1607,22 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {
                 }
             }
         }
+        if (pSPARC->cell_typ == 0) {
+            // check if axes are aligned with xyz
+            // enforce to triclinic if no aligned with xyz
+            for (j = 0; j < 3; j++) {
+                for (i = 0; i < 3; i++) {
+                    if (i == j) continue;
+                    if (fabs(pSPARC->LatVec[i+3*j]) > TEMP_TOL) {
+                        pSPARC->cell_typ = 17;
+                        i = j = 3;
+                    }
+                }
+            }
+            if (pSPARC->cell_typ > 0 && !rank) {            
+                printf(YEL"\nWARNING: This system is cuboidal. To get the best performance, please align the lattice vectors onto standard cartesian coordinate.\n" RESET);
+            }
+        }
         if (pSPARC->BC > 0) {
             if (pSPARC->BC == 1) {
                 // dirichlet boundary
@@ -1964,7 +1978,7 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {
         pSPARC->ChebDegree = Mesh2ChebDegree(h_eff);
 #ifdef DEBUG
         if (!rank && h_eff < 0.1) {
-            printf("#WARNING: for mesh less than 0.1, the default Chebyshev polynomial degree might not be enought!\n");
+            printf("WARNING: for mesh less than 0.1, the default Chebyshev polynomial degree might not be enought!\n");
         }
         if (!rank) printf("h_eff = %.2f, npl = %d\n", h_eff,pSPARC->ChebDegree);
 #endif
@@ -2222,18 +2236,6 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {
         }
     }
 
-    // scf error type, 0 - default, 1 - QE (conv_thr)
-    if (pSPARC->scf_err_type != 0 && pSPARC->scf_err_type != 1) {
-        if (!rank) printf("Cannot recognize SCF error type!\n");
-        exit(EXIT_FAILURE);
-    }
-
-    // for evaluating QE scf error, we need to perform some extra calculations
-    // e.g., an extra Poisson solve, this timer keeps track of the extra time spent
-    if (pSPARC->scf_err_type == 1) {
-        pSPARC->t_qe_extra = 0.0;
-    }
-
     // default SCF tolerance based on accuracy_level
     // we use a model curve to correlate scf tolerance and energy and force accuracy
     //     log(y) = a * log(x) + b
@@ -3389,7 +3391,7 @@ void write_output_init(SPARC_OBJ *pSPARC) {
     }
 
     fprintf(output_fp,"***************************************************************************\n");
-    fprintf(output_fp,"*                       SPARC (version Oct 31, 2023)                      *\n");
+    fprintf(output_fp,"*                       SPARC (version Jan 27, 2023)                      *\n");
     fprintf(output_fp,"*   Copyright (c) 2020 Material Physics & Mechanics Group, Georgia Tech   *\n");
     fprintf(output_fp,"*           Distributed under GNU General Public License 3 (GPL)          *\n");
     fprintf(output_fp,"*                   Start time: %s                  *\n",c_time_str);
@@ -3413,6 +3415,8 @@ void write_output_init(SPARC_OBJ *pSPARC) {
             fprintf(output_fp,"%.15f %.15f %.15f \n",pSPARC->LatVec[6],pSPARC->LatVec[7],pSPARC->LatVec[8]);
         }
     }
+    print_orthogonal_warning(pSPARC, output_fp);
+    
     if(pSPARC->cell_typ > 20 && pSPARC->cell_typ < 30)
         fprintf(output_fp,"TWIST_ANGLE: %f \n",pSPARC->twist);
     fprintf(output_fp,"FD_GRID: %d %d %d\n",pSPARC->numIntervals_x,pSPARC->numIntervals_y,pSPARC->numIntervals_z);
@@ -3462,7 +3466,7 @@ void write_output_init(SPARC_OBJ *pSPARC) {
         fprintf(output_fp,"NSTATES: %d\n",pSPARC->Nstates);
         // this should depend on temperature and preconditoner used
         if (pSPARC->Nstates < (int)(1.2*(pSPARC->Nelectron/2)+5)*pSPARC->Nspinor ) { // with kerker a factor of 1.1 might be needed
-            printf("#WARNING: Number of bands may be insufficient for efficient SCF convergence.\n");
+            printf("WARNING: Number of bands may be insufficient for efficient SCF convergence.\n");
         }
         fprintf(output_fp,"CHEB_DEGREE: %d\n",pSPARC->ChebDegree);
         if (pSPARC->CheFSI_Optmz == 1) {
@@ -3559,17 +3563,7 @@ void write_output_init(SPARC_OBJ *pSPARC) {
     fprintf(output_fp,"MAXIT_SCF: %d\n",pSPARC->MAXIT_SCF);
     fprintf(output_fp,"MINIT_SCF: %d\n",pSPARC->MINIT_SCF);
     fprintf(output_fp,"MAXIT_POISSON: %d\n",pSPARC->MAXIT_POISSON);
-    if (pSPARC->scf_err_type == 0) {
-        fprintf(output_fp,"TOL_SCF: %.2E\n",pSPARC->TOL_SCF);
-    } else if (pSPARC->scf_err_type == 1) {
-        fprintf(output_fp,"TOL_SCF_QE: %.2E\n",pSPARC->TOL_SCF);
-        if (pSPARC->spin_typ) {
-            fprintf(output_fp,"#WARNING: TOL_SCF_QE is not appropriatly set up for spin-polarized systems\n");
-        }
-        if (pSPARC->MixingVariable == 1) {
-            fprintf(output_fp,"#WARNING: TOL_SCF_QE is not appropriatly set up for potential mixing\n");
-        }
-    }
+    fprintf(output_fp,"TOL_SCF: %.2E\n",pSPARC->TOL_SCF);
     if (pSPARC->POISSON_SOLVER == 0){
         fprintf(output_fp,"POISSON_SOLVER: AAR\n");
     }else{
@@ -3615,7 +3609,7 @@ void write_output_init(SPARC_OBJ *pSPARC) {
         L_diag = sqrt(Lx * Lx + Ly * Ly + Lz * Lz);
         if (L_diag > 20.0 && pSPARC->MixingPrecond == 0) {
             fprintf(output_fp,
-                   "#WARNING: the preconditioner for SCF has been turned off, this \n"
+                   "WARNING: the preconditioner for SCF has been turned off, this \n"
                    "#might lead to slow SCF convergence. To specify SCF preconditioner, \n"
                    "#use 'MIXING_PRECOND' in the .inpt file\n");
         }
@@ -3785,6 +3779,9 @@ void write_output_init(SPARC_OBJ *pSPARC) {
             fprintf(output_fp,"EIG_PARAL_MAXNP: %d\n",pSPARC->eig_paral_maxnp);
         }
     }
+    if (pSPARC->useDefaultParalFlag == 0) {
+        fprintf(output_fp,"WARNING: Default parallelization not used. This could result in degradation of performance.\n");
+    }
     fprintf(output_fp,"***************************************************************************\n");
     fprintf(output_fp,"                             Initialization                                \n");
     fprintf(output_fp,"***************************************************************************\n");
@@ -3851,7 +3848,10 @@ void write_output_init(SPARC_OBJ *pSPARC) {
     fprintf(output_fp,"Estimated total memory usage       :  %s\n",mem_str);
     formatBytes(pSPARC->memory_usage/nproc,32,mem_str);
     fprintf(output_fp,"Estimated memory per processor     :  %s\n",mem_str);
-
+    if (pSPARC->isRbOut[0] || pSPARC->isRbOut[1] || pSPARC->isRbOut[2]) {
+        fprintf(output_fp, "WARNING: Atoms are too close to boundary for b calculation.\n");
+    }
+    
     fclose(output_fp);
 
     // write .static file
@@ -3884,6 +3884,16 @@ void write_output_init(SPARC_OBJ *pSPARC) {
 }
 
 
+void print_orthogonal_warning(SPARC_OBJ *pSPARC, FILE *output_fp) 
+{
+    printf("here\ncell_typ %d\n", pSPARC->cell_typ);
+    if (pSPARC->cell_typ == 0 || pSPARC->cell_typ >= 20) return;
+    printf("lapcT %.3e, %.3e, %.3e\n", pSPARC->lapcT[1], pSPARC->lapcT[2], pSPARC->lapcT[5]);
+    if(fabs(pSPARC->lapcT[1]) > TEMP_TOL || fabs(pSPARC->lapcT[2]) > TEMP_TOL || fabs(pSPARC->lapcT[5]) > TEMP_TOL) return;
+    printf("here\n");
+    fprintf(output_fp,"WARNING: This system is cuboidal. To get the best performance, please align the lattice vectors onto standard cartesian coordinate.\n");
+}
+
 
 /**
  * @brief Create MPI struct type SPARC_INPUT_MPI for broadcasting.
@@ -3911,7 +3921,7 @@ void SPARC_Input_MPI_create(MPI_Datatype *pSPARC_INPUT_MPI) {
                                          MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_INT, 
                                          MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_INT, 
                                          MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_INT,
-                                         MPI_INT, MPI_INT, MPI_INT, MPI_INT,
+                                         MPI_INT, MPI_INT, MPI_INT,
                                          MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
                                          MPI_DOUBLE,
                                          MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
@@ -3948,7 +3958,7 @@ void SPARC_Input_MPI_create(MPI_Datatype *pSPARC_INPUT_MPI) {
                           1, 1, 1, 1, 1, 
                           1, 1, 1, 1, 1, 
                           1, 1, 1, 1, 1,
-                          1, 1, 1, 1,/* int */ 
+                          1, 1, 1, /* int */ 
                           9, 3, L_QMASS, L_kpoint, L_kpoint,
                           L_kpoint, /* double array */
                           1, 1, 1, 1, 1, 
@@ -4012,7 +4022,6 @@ void SPARC_Input_MPI_create(MPI_Datatype *pSPARC_INPUT_MPI) {
     MPI_Get_address(&sparc_input_tmp.Nchefsi, addr + i++);
     MPI_Get_address(&sparc_input_tmp.FixRandSeed, addr + i++);
     MPI_Get_address(&sparc_input_tmp.accuracy_level, addr + i++);
-    MPI_Get_address(&sparc_input_tmp.scf_err_type, addr + i++);
     MPI_Get_address(&sparc_input_tmp.MAXIT_SCF, addr + i++);
     MPI_Get_address(&sparc_input_tmp.MINIT_SCF, addr + i++);
     MPI_Get_address(&sparc_input_tmp.MAXIT_POISSON, addr + i++);
diff --git a/src/md.c b/src/md.c
index 3871709d..35f02e53 100644
--- a/src/md.c
+++ b/src/md.c
@@ -1225,7 +1225,7 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC)
         pSPARC->ChebDegree = Mesh2ChebDegree(h_eff);
 #ifdef DEBUG
         if (!rank && h_eff < 0.1) {
-            printf("#WARNING: for mesh less than 0.1, the default Chebyshev polynomial degree might not be enought!\n");
+            printf("WARNING: for mesh less than 0.1, the default Chebyshev polynomial degree might not be enought!\n");
         }
         if (!rank) printf("h_eff = %.2f, npl = %d\n", h_eff,pSPARC->ChebDegree);
 #endif
diff --git a/src/parallelization.c b/src/parallelization.c
index f8191d09..d9334e45 100644
--- a/src/parallelization.c
+++ b/src/parallelization.c
@@ -80,6 +80,10 @@ void Setup_Comms(SPARC_OBJ *pSPARC) {
         pSPARC->npNdx = dims[0];
         pSPARC->npNdy = dims[1];
         pSPARC->npNdz = dims[2];
+        pSPARC->useDefaultParalFlag = 1;
+    } else {
+        if (!rank) printf("WARNING: Default parallelization not used. This could result in degradation of performance.\n");
+        pSPARC->useDefaultParalFlag = 0;
     }
 
     //------------------------------------------------//
@@ -1105,14 +1109,6 @@ void Setup_Comms(SPARC_OBJ *pSPARC) {
             assert(pSPARC->electronDens_in != NULL);
         }
 
-        // The following rho_in and phi_in are only used for evaluating QE scf errors
-        if (pSPARC->scf_err_type == 1) {
-            pSPARC->rho_dmcomm_phi_in = (double *)malloc(DMnd * sizeof(double));
-            assert(pSPARC->rho_dmcomm_phi_in != NULL);
-            pSPARC->phi_dmcomm_phi_in = (double *)malloc(DMnd * sizeof(double));
-            assert(pSPARC->phi_dmcomm_phi_in != NULL);
-        }
-
         // initialize electrostatic potential as random guess vector
         if (pSPARC->FixRandSeed == 1) {
             SeededRandVec(pSPARC->elecstPotential, pSPARC->DMVertices, gridsizes, -1.0, 1.0, 0);
diff --git a/src/readfiles.c b/src/readfiles.c
index 0db9166b..15657179 100644
--- a/src/readfiles.c
+++ b/src/readfiles.c
@@ -418,15 +418,9 @@ void read_input(SPARC_INPUT_OBJ *pSPARC_Input, SPARC_OBJ *pSPARC) {
             fscanf(input_fp, "%*[^\n]\n");
         } else if (strcmpi(str,"TOL_SCF:") == 0) { 
             fscanf(input_fp,"%lf",&pSPARC_Input->TOL_SCF);
-            pSPARC_Input->scf_err_type = 0; // can remove since default is 0
             snprintf(str, L_STRING, "undefined");    // initialize str
             fscanf(input_fp, "%*[^\n]\n");
-        } /*else if (strcmpi(str,"TOL_SCF_QE:") == 0) { 
-            fscanf(input_fp,"%lf",&pSPARC_Input->TOL_SCF);
-            pSPARC_Input->scf_err_type = 1;
-            snprintf(str, L_STRING, "undefined");    // initialize str
-            fscanf(input_fp, "%*[^\n]\n");
-        } */else if (strcmpi(str,"TOL_POISSON:") == 0) {
+        } else if (strcmpi(str,"TOL_POISSON:") == 0) {
             fscanf(input_fp,"%lf",&pSPARC_Input->TOL_POISSON);
             fscanf(input_fp, "%*[^\n]\n");
         } else if (strcmpi(str,"TOL_RELAX:") == 0) {
diff --git a/src/relax.c b/src/relax.c
index 321d0ff8..d1770b91 100644
--- a/src/relax.c
+++ b/src/relax.c
@@ -1236,7 +1236,7 @@ void reinitialize_cell_mesh(SPARC_OBJ *pSPARC, double vol)
             pSPARC->ChebDegree = Mesh2ChebDegree(h_eff);
     #ifdef DEBUG
             if (!rank && h_eff < 0.1) {
-                printf("#WARNING: for mesh less than 0.1, the default Chebyshev polynomial degree might not be enought!\n");
+                printf("WARNING: for mesh less than 0.1, the default Chebyshev polynomial degree might not be enought!\n");
             }
             if (!rank) printf("h_eff = %.2f, npl = %d\n", h_eff,pSPARC->ChebDegree);
     #endif