Merge pull request #196 from YaphetS-jx/dev

ChangeLog

@@ -3,6 +3,14 @@
 -Name
 -changes
 
+--------------
+Sep 18, 2023
+Name: Xin Jing
+Changes: (highT/)
+1. Distinguishing between SPARC and dev_SPARC
+2. Change a few spin tests to include NP_SPIN_PARAL: 1
+3. Fix a bug in initial electron density and magnetization
+
 --------------
 Sep 08, 2023
 Name: Xin Jing

src/eigenSolver.c

@@ -354,9 +354,12 @@ void CheFSI(SPARC_OBJ *pSPARC, double lambda_cutoff, double *x0, int count, int
     #else
     // allocate memory for block cyclic format of the wavefunction
     if (pSPARC->npband > 1 || pSPARC->Nspinor_eig != pSPARC->Nspinor_spincomm) {
-        pSPARC->Yorb_BLCYC = (double *)malloc(
-            pSPARC->nr_orb_BLCYC * pSPARC->nc_orb_BLCYC * sizeof(double));
-        assert(pSPARC->Yorb_BLCYC != NULL);
+        pSPARC->Xorb_BLCYC = (double *)malloc(pSPARC->nr_orb_BLCYC * pSPARC->nc_orb_BLCYC * sizeof(double));
+        pSPARC->Yorb_BLCYC = (double *)malloc(pSPARC->nr_orb_BLCYC * pSPARC->nc_orb_BLCYC * sizeof(double));
+        assert(pSPARC->Xorb_BLCYC != NULL && pSPARC->Yorb_BLCYC != NULL);
+    } else {
+        pSPARC->Xorb_BLCYC = pSPARC->Xorb + spn_i*DMnd;
+        pSPARC->Yorb_BLCYC = pSPARC->Yorb + spn_i*DMnd;
     }
     Project_Hamiltonian(pSPARC, pSPARC->DMVertices_dmcomm, pSPARC->Yorb + spn_i*DMnd, DMndsp, pSPARC->Xorb + spn_i*DMnd, DMndsp,
                         pSPARC->Hp, pSPARC->Mp, k, spn_i, pSPARC->dmcomm);
@@ -416,14 +419,6 @@ void CheFSI(SPARC_OBJ *pSPARC, double lambda_cutoff, double *x0, int count, int
     #ifdef USE_DP_SUBEIG
     DP_Subspace_Rotation(pSPARC, pSPARC->Xorb + spn_i*DMnd);
     #else
-    if (pSPARC->npband > 1 || pSPARC->Nspinor_eig != pSPARC->Nspinor_spincomm) {
-        // find Y * Q, store the result in Xorb (band+domain) and Xorb_BLCYC (block cyclic format)
-        pSPARC->Xorb_BLCYC = (double *)malloc(pSPARC->nr_orb_BLCYC * pSPARC->nc_orb_BLCYC * sizeof(double));
-        assert(pSPARC->Xorb_BLCYC != NULL);
-    } else {
-        pSPARC->Xorb_BLCYC = pSPARC->Xorb + spn_i*DMnd;
-    }
-
     // find Y * Q, store the result in Xorb (band+domain) and Xorb_BLCYC (block cyclic format)
     Subspace_Rotation(pSPARC, pSPARC->Yorb_BLCYC, pSPARC->Q, 
                         pSPARC->Xorb_BLCYC, pSPARC->Xorb + spn_i*DMnd, k, spn_i);
@@ -1505,9 +1500,7 @@ void Project_Hamiltonian(SPARC_OBJ *pSPARC, int *DMVertices, double *Y, int ldi,
         // distribute orbitals into block cyclic format
         pdgemr2d_(&DMndspe, &pSPARC->Nstates, Y, &ONE, &ONE, pSPARC->desc_orbitals,
                   pSPARC->Yorb_BLCYC, &ONE, &ONE, pSPARC->desc_orb_BLCYC, &pSPARC->ictxt_blacs); 
-    } else {
-        pSPARC->Yorb_BLCYC = Y;
-    }
+    } 
     t2 = MPI_Wtime();  
     #ifdef DEBUG  
     if(!rank && spn_i == 0) 
@@ -1528,7 +1521,7 @@ void Project_Hamiltonian(SPARC_OBJ *pSPARC, int *DMVertices, double *Y, int ldi,
                 &ONE, &ONE, pSPARC->desc_Mp_BLCYC);
         } else {
             #ifdef DEBUG    
-            if (!rank && spn_i == 0) printf("rank = %d, STARTING PDGEMM ...\n",rank);
+            if (!rank && spn_i == 0) printf("rank = %d, STARTING PDSYRK ...\n",rank);
             #endif   
             // perform matrix multiplication using ScaLAPACK routines
             pdsyrk_("U", "T", &pSPARC->Nstates, &DMndspe, &alpha, pSPARC->Yorb_BLCYC, &ONE, &ONE,

src/electronicGroundState.c

@@ -84,13 +84,13 @@ void Calculate_electronicGroundState(SPARC_OBJ *pSPARC) {
         else
             SCF_ind = 1;
         if (pSPARC->REFERENCE_CUTOFF > 0.5*nn) {
-            printf("\nWARNING: REFERENCE _CUFOFF (%.6f Bohr) > 1/2 nn (nearest neighbor) distance (%.6f Bohr) in SCF#%d\n",
+            printf("\nWARNING: REFERENCE_CUFOFF (%.6f Bohr) > 1/2 nn (nearest neighbor) distance (%.6f Bohr) in SCF#%d\n",
                         pSPARC->REFERENCE_CUTOFF, 0.5*nn,  SCF_ind);
         }
         if (pSPARC->REFERENCE_CUTOFF < pSPARC->delta_x ||
             pSPARC->REFERENCE_CUTOFF < pSPARC->delta_y ||
             pSPARC->REFERENCE_CUTOFF < pSPARC->delta_z ) {
-            printf("\nWARNING: REFERENCE _CUFOFF (%.6f Bohr) < MESH_SPACING (dx %.6f Bohr, dy %.6f Bohr, dz %.6f Bohr) in SCF#%d\n",
+            printf("\nWARNING: REFERENCE_CUFOFF (%.6f Bohr) < MESH_SPACING (dx %.6f Bohr, dy %.6f Bohr, dz %.6f Bohr) in SCF#%d\n",
                         pSPARC->REFERENCE_CUTOFF, pSPARC->delta_x, pSPARC->delta_y, pSPARC->delta_z, SCF_ind );
         }
     }
@@ -247,8 +247,7 @@ void Calculate_electronicGroundState(SPARC_OBJ *pSPARC) {
         }
     } else if(pSPARC->Calc_pres == 1){
         t1 = MPI_Wtime();
-        Calculate_electronic_pressure(pSPARC);
-        // if (pSPARC->d3Flag == 1) d3_grad_cell_stress(pSPARC); // add the D3 contribution on pressure to total pressure. move this function into Calculate_electronic_pressure?
+        Calculate_electronic_pressure(pSPARC);        
         t2 = MPI_Wtime();
         if(!rank && pSPARC->Verbosity) {
         	output_fp = fopen(pSPARC->OutFilename,"a");
@@ -583,7 +582,7 @@ void scf(SPARC_OBJ *pSPARC)
             Transfer_Veff_loc(pSPARC, pSPARC->Veff_loc_dmcomm_phi + i*DMnd, pSPARC->Veff_loc_dmcomm + i*pSPARC->Nd_d_dmcomm);
     }
 
-	#ifdef DEBUG
+    #ifdef DEBUG
     t2 = MPI_Wtime();
     if (rank == 0) {
         printf("rank = %d, Veff calculation and Bcast (non-blocking) took %.3f ms\n",rank,(t2-t1)*1e3); 
@@ -608,13 +607,6 @@ void scf_loop(SPARC_OBJ *pSPARC) {
 
     int DMnd = pSPARC->Nd_d;    
     int i, k, SCFcount;
-
-#ifdef DEBUG
-    int spn_i;
-    int Nk = pSPARC->Nkpts_kptcomm;
-    int Ns = pSPARC->Nstates;
-#endif
-
     double error, dEtot, dEband;
     double t_scf_s, t_scf_e, t_cum_scf;
     double Veff_mean[4];
@@ -914,6 +906,9 @@ void scf_loop(SPARC_OBJ *pSPARC) {
         }
     } else {
         #ifdef DEBUG
+        int spn_i;
+        int Nk = pSPARC->Nkpts_kptcomm;
+        int Ns = pSPARC->Nstates;
         int spin_maxocc = 0, k_maxocc = 0;
         double maxocc = -1.0;
         for (spn_i = 0; spn_i < pSPARC->Nspin_spincomm; spn_i++){
@@ -928,7 +923,6 @@ void scf_loop(SPARC_OBJ *pSPARC) {
                     k_maxocc = k;
                 }
 
-                // #ifdef DEBUG
                 if(!rank) {
                     int nocc_print = min(200,pSPARC->Nstates - pSPARC->Nelectron/2 + 10);
                     nocc_print = min(nocc_print, pSPARC->Nstates);
@@ -964,8 +958,8 @@ void scf_loop(SPARC_OBJ *pSPARC) {
             if (pSPARC->BC != 1) printf("\nk = [%.3f, %.3f, %.3f]\n", k1_red, k2_red, k3_red);
             printf("Occupation of state %d (90%%) = %.15f.\n"
                 "Occupation of state %d (100%%) = %.15f.\n",
-                ind_90percent+1, (3.0-pSPARC->Nspin)/pSPARC->Nspinor * g_ind_90percent,
-                ind_100percent+1, (3.0-pSPARC->Nspin)/pSPARC->Nspinor * g_ind_100percent);
+                ind_90percent+1, pSPARC->occfac * g_ind_90percent,
+                ind_100percent+1, pSPARC->occfac * g_ind_100percent);
         }
         #endif
     }

src/electrostatics.c

@@ -47,7 +47,7 @@ void Calculate_PseudochargeCutoff(SPARC_OBJ *pSPARC) {
 #define V_temp(i,j,k) V_temp[(k)*(xlen_ex)*(ylen_ex)+(j)*(xlen_ex)+(i)]
 // TODO: rb_max depends on finite difference order as well, improve this function
 #define RB_MAX(h) (h) < 1.5 ? (h)*10+10 : 20*(h)-9.5
-#define RB_Y(y,L,typ) (typ) < 20 ? (y) : acos(1 - (y)*(y)/(2*(L)*(L))) + pSPARC->twist*y
+#define RB_Y(y,L,typ) (typ) < 20 ? (y) : asin(y/L) + 0.5*pSPARC->twist*y //acos(1 - (y)*(y)/(2*(L)*(L))) + pSPARC->twist*y
 //#define RB_MAX(h) (h) < 1.5 ? (h)*30+5.5 : 40*(h)-9.5 
     double Rbmax_x, Rbmax_y, Rbmax_z;
     int rank, nx, ny, nz, nxp, nyp, nzp, FDn, rlen_ex, rlen, i, j, k, ityp, xlen, ylen, zlen, xlen_ex, ylen_ex, zlen_ex;
@@ -769,7 +769,7 @@ void GetInfluencingAtoms(SPARC_OBJ *pSPARC) {
  * @brief   Calculate pseudocharge density for given atom positions
  */
 void Generate_PseudoChargeDensity(SPARC_OBJ *pSPARC) {
-#define electronDens_at(s,i,j,k) electronDens_at[s+(k)*DMnx*DMny+(j)*DMnx+(i)]
+#define electronDens_at(i,j,k) electronDens_at[(k)*DMnx*DMny+(j)*DMnx+(i)]
 #define electronDens_core(i,j,k) electronDens_core[(k)*DMnx*DMny+(j)*DMnx+(i)]
 #define psdChrgDens(i,j,k) psdChrgDens[(k)*DMnx*DMny+(j)*DMnx+(i)]
 #define psdChrgDens_ref(i,j,k) psdChrgDens_ref[(k)*DMnx*DMny+(j)*DMnx+(i)]
@@ -839,7 +839,7 @@ void Generate_PseudoChargeDensity(SPARC_OBJ *pSPARC) {
     DMnd = pSPARC->Nd_d;
 
     // initialize to zero at the beginning of each relax/MD step
-    memset(pSPARC->electronDens_at, 0, sizeof(double)*DMnd*pSPARC->Nspdentd);
+    memset(pSPARC->electronDens_at, 0, sizeof(double)*DMnd);
     memset(pSPARC->electronDens_core, 0, sizeof(double)*DMnd);
     memset(pSPARC->psdChrgDens, 0, sizeof(double)*DMnd);
     memset(pSPARC->psdChrgDens_ref, 0, sizeof(double)*DMnd);
@@ -848,13 +848,13 @@ void Generate_PseudoChargeDensity(SPARC_OBJ *pSPARC) {
     double *magx, *magy, *magz;
     magx = magy = magz = NULL;
     if (pSPARC->spin_typ == 1) {
-        magz = pSPARC->mag;
-        memset(pSPARC->mag, 0, sizeof(double)*DMnd*pSPARC->Nmag);
+        memset(pSPARC->mag_at, 0, sizeof(double)*DMnd);
+        magz = pSPARC->mag_at;
     } else if (pSPARC->spin_typ == 2) {
-        magx = pSPARC->mag + DMnd;
-        magy = pSPARC->mag + DMnd*2;
-        magz = pSPARC->mag + DMnd*3;
-        memset(pSPARC->mag, 0, sizeof(double)*DMnd*pSPARC->Nmag);
+        memset(pSPARC->mag_at, 0, sizeof(double)*DMnd*3);
+        magx = pSPARC->mag_at;
+        magy = pSPARC->mag_at + DMnd;
+        magz = pSPARC->mag_at + DMnd*2;
     }
 
     // calculate pseudocharge density bJ and (self + correction) energy
@@ -1027,9 +1027,9 @@ void Generate_PseudoChargeDensity(SPARC_OBJ *pSPARC) {
                         i_global = i + pSPARC->Atom_Influence_local[ityp].xs[iat];
                         i_DM = i_global - pSPARC->DMVertices[0]; // local coord 
                         if (i_DM < 0 || i_DM >= DMnx) continue;
-                        pSPARC->electronDens_at(0,i_DM,j_DM,k_DM) += rho_J(ip,jp,kp);
+                        pSPARC->electronDens_at(i_DM,j_DM,k_DM) += rho_J(ip,jp,kp);
                         pSPARC->electronDens_core(i_DM,j_DM,k_DM) += rho_c_J(ip,jp,kp);
-                        if (pSPARC->spin_typ == 1) {                            
+                        if (pSPARC->spin_typ == 1) {
                             magz(i_DM,j_DM,k_DM) += (pSPARC->Atom_Influence_local[ityp].atom_spin[3*iat+2] / pSPARC->Znucl[ityp]) * rho_J(ip,jp,kp);
                         } else if (pSPARC->spin_typ == 2) {
                             magx(i_DM,j_DM,k_DM) += (pSPARC->Atom_Influence_local[ityp].atom_spin[3*iat] / pSPARC->Znucl[ityp]) * rho_J(ip,jp,kp);
@@ -1093,13 +1093,6 @@ void Generate_PseudoChargeDensity(SPARC_OBJ *pSPARC) {
         }
     }
 
-    if (pSPARC->spin_typ == 1) {
-        Calculate_diagonal_Density(pSPARC, DMnd, magz, pSPARC->electronDens_at, pSPARC->electronDens_at+DMnd, pSPARC->electronDens_at+2*DMnd);
-    } else if (pSPARC->spin_typ == 2) {
-        Calculate_Magnorm(pSPARC, DMnd, magx, magy, magz, pSPARC->mag);
-        Calculate_diagonal_Density(pSPARC, DMnd, pSPARC->mag, pSPARC->electronDens_at, pSPARC->electronDens_at+DMnd, pSPARC->electronDens_at+2*DMnd);
-    }    
-
     /*  Calculate integral of b and Esc  */
     double int_b = 0.0, int_rho = 0.0;
     // find integral of b, Esc locally
@@ -1109,7 +1102,7 @@ void Generate_PseudoChargeDensity(SPARC_OBJ *pSPARC) {
             for (j = 0; j < DMny; j++) {
                 for (i = 0; i < DMnx; i++) {
                     int_b += pSPARC->psdChrgDens(i,j,k) * pSPARC->Intgwt_phi[count];
-                    int_rho += pSPARC->electronDens_at(0,i,j,k) * pSPARC->Intgwt_phi[count];
+                    int_rho += pSPARC->electronDens_at(i,j,k) * pSPARC->Intgwt_phi[count];
                     Esc += (pSPARC->psdChrgDens(i,j,k) + pSPARC->psdChrgDens_ref(i,j,k)) * pSPARC->Vc(i,j,k) * pSPARC->Intgwt_phi[count];
                     count++;
                 }
@@ -1121,7 +1114,7 @@ void Generate_PseudoChargeDensity(SPARC_OBJ *pSPARC) {
             for (j = 0; j < DMny; j++) {
                 for (i = 0; i < DMnx; i++) {
                     int_b += pSPARC->psdChrgDens(i,j,k);
-                    int_rho += pSPARC->electronDens_at(0,i,j,k);
+                    int_rho += pSPARC->electronDens_at(i,j,k);
                     Esc += (pSPARC->psdChrgDens(i,j,k) + pSPARC->psdChrgDens_ref(i,j,k)) * pSPARC->Vc(i,j,k);
                 }
             }
@@ -1153,16 +1146,12 @@ void Generate_PseudoChargeDensity(SPARC_OBJ *pSPARC) {
     if (pSPARC->CyclixFlag) {
         for (int i = 0; i < DMnd; i++) {
             pSPARC->electronDens_at[i] *= scal_fac;
-            for(int n = 1; n < pSPARC->Nspdentd; n++) 
-                pSPARC->electronDens_at[n*DMnd + i] *= scal_fac;
             Nelectron_check += pSPARC->electronDens_at[i] * pSPARC->Intgwt_phi[i];            
         }                    
     } else {
         for (int i = 0; i < DMnd; i++) {
             pSPARC->electronDens_at[i] *= scal_fac;
             Nelectron_check += pSPARC->electronDens_at[i];
-            for(int n = 1; n < pSPARC->Nspdentd; n++) 
-                pSPARC->electronDens_at[n*DMnd + i] *= scal_fac;
         }
         Nelectron_check *= pSPARC->dV;
     }

src/finalization.c

@@ -169,6 +169,7 @@ void Free_basic(SPARC_OBJ *pSPARC) {
         free(pSPARC->electronDens);
         if (pSPARC->spin_typ > 0) {
             free(pSPARC->mag);
+            free(pSPARC->mag_at);
             free(pSPARC->AtomMag);
         }
         free(pSPARC->psdChrgDens);
@@ -181,6 +182,26 @@ void Free_basic(SPARC_OBJ *pSPARC) {
         }    
         free(pSPARC->elecstPotential);
         free(pSPARC->Veff_loc_dmcomm_phi);
+
+        // free MD and relax stuff
+    	if(pSPARC->MDFlag == 1 || pSPARC->RelaxFlag == 1 || pSPARC->RelaxFlag == 3){
+      		free(pSPARC->delectronDens);
+      		free(pSPARC->delectronDens_0dt);
+      		free(pSPARC->delectronDens_1dt);
+      		free(pSPARC->delectronDens_2dt);
+      		free(pSPARC->atom_pos_nm);
+      		free(pSPARC->atom_pos_0dt);
+      		free(pSPARC->atom_pos_1dt);
+      		free(pSPARC->atom_pos_2dt);
+    	}
+
+        if (pSPARC->spin_typ == 2) {
+            free(pSPARC->XCPotential_nc);
+        }
+    }
+
+    // mixing variables
+    if (pSPARC->dmcomm_phi != MPI_COMM_NULL) {
         free(pSPARC->mixing_hist_xk);
         free(pSPARC->mixing_hist_fk);
         free(pSPARC->mixing_hist_fkm1);
@@ -191,6 +212,9 @@ void Free_basic(SPARC_OBJ *pSPARC) {
         // saving potential history
         if (pSPARC->MixingVariable == 1) {
             free(pSPARC->Veff_loc_dmcomm_phi_in);
+            if (pSPARC->spin_typ == 2) {
+                free(pSPARC->Veff_dia_loc_dmcomm_phi);
+            }
         }
 
         if (pSPARC->MixingVariable == 0 && pSPARC->spin_typ) {
@@ -204,25 +228,6 @@ void Free_basic(SPARC_OBJ *pSPARC) {
         }
 
         free(pSPARC->mixing_hist_Pfk);
-
-        // free MD and relax stuff
-    	if(pSPARC->MDFlag == 1 || pSPARC->RelaxFlag == 1 || pSPARC->RelaxFlag == 3){
-      		free(pSPARC->delectronDens);
-      		free(pSPARC->delectronDens_0dt);
-      		free(pSPARC->delectronDens_1dt);
-      		free(pSPARC->delectronDens_2dt);
-      		free(pSPARC->atom_pos_nm);
-      		free(pSPARC->atom_pos_0dt);
-      		free(pSPARC->atom_pos_1dt);
-      		free(pSPARC->atom_pos_2dt);
-    	}
-
-        if (pSPARC->spin_typ == 2) {
-            free(pSPARC->XCPotential_nc);
-            if (pSPARC->MixingVariable == 1) {
-                free(pSPARC->Veff_dia_loc_dmcomm_phi);
-            }
-        }
     }
 
     // free preconditioner coeff arrays
@@ -280,18 +285,6 @@ void Free_SPARC(SPARC_OBJ *pSPARC) {
         free(pSPARC->IP_displ_SOC); 
 
     if (pSPARC->usefock > 0) {
-        free(pSPARC->k1_hf);
-        free(pSPARC->k2_hf);
-        free(pSPARC->k3_hf);
-        free(pSPARC->kpthf_ind);
-        free(pSPARC->kpthf_ind_red);
-        free(pSPARC->kpthfred2kpthf);
-        free(pSPARC->kpthf_pn);
-        free(pSPARC->kpts_hf_red_list);
-        free(pSPARC->k1_shift);
-        free(pSPARC->k2_shift);
-        free(pSPARC->k3_shift);
-        free(pSPARC->Kptshift_map);
         free_exx(pSPARC);
     }
 
@@ -332,7 +325,7 @@ void Free_SPARC(SPARC_OBJ *pSPARC) {
     if (pSPARC->ixc[3] != 0){
         vdWDF_free(pSPARC);
     }
-    if(pSPARC->ixc[2] == 1) {
+    if (pSPARC->ixc[2] == 1) {
         free_MGGA(pSPARC);
     }
 

src/include/eigenSolver.h

@@ -236,6 +236,13 @@ int Mesh2ChebDegree(double h);
 
 /**
  * @brief   Orthogonalization of dense matrix A by Choleskey factorization
+ * 
+ * @param A            (INPUT)  Distributed dense matrix A.
+ * @param descA        (INPUT)  Descriptor of A.
+ * @param z            (INPUT/OUTPUT) INPUT: z=A'*A, OUTPUT: A'*A=z'*z, z is upper triangular matrix.
+ * @param descz        (INPUT)  Descriptor of Z.
+ * @param m            (INPUT)  Row blocking factor.
+ * @param n            (INPUT)  Column blocking factor.
  */
 void Chol_orth(double *A, const int *descA, double *z, const int *descz, const int *m, const int *n);
 

src/include/forces.h

@@ -20,7 +20,6 @@
 /**
  * @brief    Calculate atomic forces.
  */
-//void Calculate_Atomicforces(SPARC_OBJ *pSPARC);
 void Calculate_EGS_Forces(SPARC_OBJ *pSPARC);
 
 /**