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.

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