Firx bug in non-orthogonal cell and add warnings

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

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
 **/

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");
+        }
     }
 
 }

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);
     }

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);

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.
  *

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