Merge pull request SPARC-X#190 from StevenZhangCSFM/add_rSCAN_r2SCAN

ChangeLog

@@ -3,6 +3,14 @@
 -Name
 -changes
 
+--------------
+Aug 01, 2023
+Name: Boqin Zhang
+Changes: (exchangeCorrelation.c, xc/mgga/, md.c)
+1. Add rscan metaGGA functional
+2. Remove vdWDFflag and mGGAflag. Replaced by ixc[2], ixc[3]
+3. Fix the bug in calculating minimum distance in md.c
+
 --------------
 Jul 10, 2023
 Name: Boqin Zhang

doc/.LaTeX/System.tex

@@ -329,13 +329,13 @@
 Choice of exchange-correlation functional. Options are \texttt{LDA\_PW} (Perdew-Wang LDA), \texttt{LDA\_PZ} (Purdew-Zunger LDA), 
 \texttt{GGA\_PBE} (PBE GGA), \texttt{GGA\_RPBE} (revised PBE GGA), and \texttt{GGA\_PBEsol} (PBE GGA revised for solids), 
 \texttt{PBE0}, \texttt{HF} (Hartree-Fock), \texttt{HSE},
-\texttt{vdWDF1} (van der Waals Density Functional developed by Dion et al.), \texttt{vdWDF2} (vdW Density Functional modified by Lee et al), and \texttt{SCAN} (SCAN metaGGA), and \texttt{R2SCAN} (r2SCAN metaGGA).
+\texttt{vdWDF1} (van der Waals Density Functional developed by Dion et al.), \texttt{vdWDF2} (vdW Density Functional modified by Lee et al), \texttt{SCAN} (SCAN metaGGA), \texttt{RSCAN} (rSCAN metaGGA), and \texttt{R2SCAN} (r2SCAN metaGGA).
 \end{block}
 
 \begin{block}{Remark}
 For spin-polarized calculation (\hyperlink{SPIN_TYP}{\texttt{SPIN\_TYP}} = 1), \texttt{LDA\_PZ} is not available.
 
-Currently SCAN and R2SCAN does not support nonlinear core correction pseudopotential.
+Currently SCAN, RSCAN and R2SCAN does not support nonlinear core correction pseudopotential.
 \end{block}
 
 \end{frame}

src/electronicGroundState.c

@@ -123,7 +123,7 @@ void Calculate_electronicGroundState(SPARC_OBJ *pSPARC) {
         if (pSPARC->d3Flag == 1) {
         	fprintf(output_fp,"DFT-D3 correction                  :%18.10E (Ha)\n", pSPARC->d3Energy[0]);
         }
-        if (pSPARC->vdWDFFlag != 0) {
+        if (pSPARC->ixc[3] != 0) {
             fprintf(output_fp,"vdWDF energy                       :%18.10E (Ha)\n", pSPARC->vdWDFenergy);
         }
         fclose(output_fp);

src/energy.c

@@ -90,7 +90,7 @@ void Calculate_Free_Energy(SPARC_OBJ *pSPARC, double *electronDens)
             E2 *= 0.5 * pSPARC->dV;
             E3 *= pSPARC->dV;
 
-            if (pSPARC->mGGAflag == 1) {
+            if (pSPARC->ixc[2]) {
                 double Emgga;
                 if (pSPARC->spin_typ == 0)
                     VectorDotProduct(pSPARC->KineticTauPhiDomain, pSPARC->vxcMGGA3, pSPARC->Nd_d, &Emgga, pSPARC->dmcomm_phi);

src/exchangeCorrelation.c

@@ -24,6 +24,7 @@
 #include "vdWDFnonlinearCorre.h"
 #include "mGGAtauTransferTauVxc.h"
 #include "mGGAscan.h"
+#include "mGGArscan.h"
 #include "mGGAr2scan.h"
 
 /**
@@ -99,6 +100,9 @@ void Calculate_Vxc(SPARC_OBJ *pSPARC)
         case 4:
             scanx(DMnd, rho, sigma, tau, ex, vx, v2x, v3x);
             break;
+        case 5:
+            rscanx(DMnd, rho, sigma, tau, ex, vx, v2x, v3x);
+            break;
         case 6:
             r2scanx(DMnd, rho, sigma, tau, ex, vx, v2x, v3x);
             break;
@@ -131,6 +135,9 @@ void Calculate_Vxc(SPARC_OBJ *pSPARC)
         case 4:
             scanc(DMnd, rho, sigma, tau, ec, vc, v2c, v3c);
             break;
+        case 5:
+            rscanc(DMnd, rho, sigma, tau, ec, vc, v2c, v3c);
+            break;
         case 6:
             r2scanc(DMnd, rho, sigma, tau, ec, vc, v2c, v3c);
             break;
@@ -175,7 +182,7 @@ void Calculate_Vxc(SPARC_OBJ *pSPARC)
             }
         }
 
-        if (pSPARC->ixc[3])
+        if (pSPARC->ixc[3] != 0)
             Calculate_nonLinearCorr_E_V_vdWDF(pSPARC, rho);
 
         if (pSPARC->isgradient) {
@@ -255,6 +262,9 @@ void Calculate_Vxc(SPARC_OBJ *pSPARC)
         case 4:
             scanx_spin(DMnd, rho, sigma, tau, ex, vx, v2x, v3x);
             break;
+        case 5:
+            rscanx_spin(DMnd, rho, sigma, tau, ex, vx, v2x, v3x);
+            break;
         case 6:
             r2scanx_spin(DMnd, rho, sigma, tau, ex, vx, v2x, v3x);
             break;
@@ -287,6 +297,9 @@ void Calculate_Vxc(SPARC_OBJ *pSPARC)
         case 4:
             scanc_spin(DMnd, rho, sigma, tau, ec, vc, v2c, v3c);
             break;
+        case 5:
+            rscanc_spin(DMnd, rho, sigma, tau, ec, vc, v2c, v3c);
+            break;
         case 6:
             r2scanc_spin(DMnd, rho, sigma, tau, ec, vc, v2c, v3c);
             break;
@@ -398,7 +411,12 @@ void Calculate_Vxc(SPARC_OBJ *pSPARC)
                 pSPARC->ixc[2] = 1; pSPARC->ixc[3] = 0;
                 pSPARC->xcoption[0] = 0; pSPARC->xcoption[1] = 0;
             }
-            else if (strcmpi(pSPARC->XC, "R2SCAN") == 0) {
+            else if (strcmpi(pSPARC->XC, "RSCAN") == 0) {
+                pSPARC->ixc[0] = 5; pSPARC->ixc[1] = 5; 
+                pSPARC->ixc[2] = 1; pSPARC->ixc[3] = 0;
+                pSPARC->xcoption[0] = 0; pSPARC->xcoption[1] = 0;
+            }
+            else {
                 pSPARC->ixc[0] = 6; pSPARC->ixc[1] = 6; 
                 pSPARC->ixc[2] = 1; pSPARC->ixc[3] = 0;
                 pSPARC->xcoption[0] = 0; pSPARC->xcoption[1] = 0;

src/finalization.c

@@ -313,10 +313,10 @@ void Free_SPARC(SPARC_OBJ *pSPARC) {
     if (pSPARC->d3Flag == 1) {
         free_D3_coefficients(pSPARC); // this function is moved from electronicGroundState.c
     }
-    if (pSPARC->vdWDFFlag != 0){
+    if (pSPARC->ixc[3] != 0){
         vdWDF_free(pSPARC);
     }
-    if(pSPARC->mGGAflag == 1) {
+    if(pSPARC->ixc[2] == 1) {
         free_MGGA(pSPARC);
     }
 

src/hamiltonianVecRoutines.c

@@ -90,7 +90,7 @@ void Hamiltonian_vectors_mult(
     }
 
     // adding metaGGA term
-    if(pSPARC->mGGAflag == 1 && pSPARC->countPotentialCalculate > 1) {
+    if(pSPARC->ixc[2] && pSPARC->countPotentialCalculate > 1) {
 
         int Lanczos_flag = (comm == pSPARC->kptcomm_topo) ? 1 : 0;
         int sg = pSPARC->spin_start_indx + spin;
@@ -192,7 +192,7 @@ void Hamiltonian_vectors_mult_kpt(
         #endif
 
         // adding metaGGA term
-        if(pSPARC->mGGAflag == 1 && pSPARC->countPotentialCalculate > 1) {
+        if(pSPARC->ixc[2] && pSPARC->countPotentialCalculate > 1) {
 
             int Lanczos_flag = (comm == pSPARC->kptcomm_topo) ? 1 : 0;
             int sg = pSPARC->spin_start_indx + spin;

src/include/isddft.h

@@ -764,7 +764,6 @@ typedef struct _SPARC_OBJ{
     double d3Sigma[9]; // stress of cell caused by d3
 
     // vdW-DF1 and vdW-DF2
-    int vdWDFFlag;
     int vdWDFnrpoints;
     int vdWDFnqs;
     double *vdWDFecLinear;
@@ -805,7 +804,6 @@ typedef struct _SPARC_OBJ{
     int countPotentialCalculate; // for helping output variables in 1st step, to be deleted in the future
 
     /* metaGGA functionals (SCAN) */
-    int mGGAflag;
     double *KineticTauPhiDomain;
     double *vxcMGGA1; // d(n\epsilon)/dn, in dmcomm_phi
     double *vxcMGGA2; // d(n\epsilon)/d|grad n|, in dmcomm_phi

src/initialization.c

@@ -337,13 +337,13 @@ void Initialize(SPARC_OBJ *pSPARC, int argc, char *argv[]) {
     }
 
     // initialize vdW-DF
-    if (pSPARC->vdWDFFlag != 0) {
+    if (pSPARC->ixc[3] != 0) {
         vdWDF_initial_read_kernel(pSPARC); // read kernel function and 2nd derivative of spline functions
         // printf("rank %d, d2 of kernel function vdWDFd2Phidk2[2][4]=%.9e\n", rank, pSPARC->vdWDFd2Phidk2[2][4]); // to verify it
     }
 
     // initialize metaGGA
-    if(pSPARC->mGGAflag == 1) {
+    if(pSPARC->ixc[2]) {
         initialize_MGGA(pSPARC);
     }
 
@@ -1246,21 +1246,11 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {
     pSPARC->NLCC_flag = NLCC_flag;
 
     // check if exchange-correlation functional is metaGGA
-    pSPARC->mGGAflag = 0;
-    if ((strcmpi(pSPARC->XC, "SCAN") == 0) || (strcmpi(pSPARC->XC, "R2SCAN") == 0)) { // it can be expand, such as adding r2SCAN 
+    if (pSPARC->ixc[2]) { // it can be expand, such as adding r2SCAN 
         if (pSPARC->NLCC_flag) {
-            if (!rank) printf("\nERROR: currently SCAN and R2SCAN functional does not support applying NLCC pseudopotential!\n");
+            if (!rank) printf("\nERROR: currently SCAN, RSCAN and R2SCAN functional does not support applying NLCC pseudopotential!\n");
             exit(EXIT_FAILURE);
         }
-        pSPARC->mGGAflag = 1;
-    }
-    // check if exchange-correlation functional is vdW-DF1 or vdW-DF2
-    pSPARC->vdWDFFlag = 0;
-    if (strcmpi(pSPARC->XC,"vdWDF1") == 0) {
-        pSPARC->vdWDFFlag = 1;
-    }
-    if (strcmpi(pSPARC->XC,"vdWDF2") == 0) {
-        pSPARC->vdWDFFlag = 2;
     }
 
     // initialize energy values
@@ -2220,7 +2210,7 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {
         && fabs(pSPARC->k3[0]) < TEMP_TOL
         && pSPARC->SOC_Flag == 0);
 
-    if (pSPARC->vdWDFFlag != 0){
+    if (pSPARC->ixc[3] != 0){
         if ((pSPARC->BCx)||(pSPARC->BCy)||(pSPARC->BCz)) {
             if (rank == 0)
                 printf(RED "ERROR: vdW-DF does not support Dirichlet boundary condition!\n" RESET);
@@ -2229,7 +2219,7 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {
     }
 
     #if !defined(USE_MKL) && !defined(USE_FFTW)
-    if (pSPARC->vdWDFFlag != 0){
+    if (pSPARC->ixc[3] != 0){
         if (rank == 0)
             printf(RED "ERROR: To use vdW-DF, please turn on MKL or FFTW in makefile!\n"
             "Or you can stop using vdW-DF by setting other exchange-correlation functionals.\n" RESET);
@@ -2242,7 +2232,7 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {
     }
     #endif // #if !defined(USE_MKL) && !defined(USE_FFTW)
 
-    if (pSPARC->mGGAflag == 1) {
+    if (pSPARC->ixc[2]) {
         // if (pSPARC->spin_typ != 0) {
         //     if (rank == 0) 
         //         printf(RED "ERROR: currently SCAN does not support spin polarization!\n" RESET);
@@ -2336,7 +2326,7 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {
 
     // constraints on SOC 
     if (pSPARC->SOC_Flag == 1) {
-        if (pSPARC->usefock || pSPARC->mGGAflag || pSPARC->SQFlag) {
+        if (pSPARC->usefock || pSPARC->ixc[2] || pSPARC->SQFlag) {
             if (rank == 0) 
                 printf(RED "ERROR: Hybrid functional, SCAN and SQ are not supported in this version of spin-orbit coupling implementation.\n" RESET);
             exit(EXIT_FAILURE);
@@ -2349,7 +2339,7 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {
     }
 
     if (pSPARC->CyclixFlag == 1) {
-        if (pSPARC->usefock || pSPARC->SQFlag || pSPARC->mGGAflag) {
+        if (pSPARC->usefock || pSPARC->SQFlag || pSPARC->ixc[2]) {
             if (rank == 0) 
                 printf(RED "ERROR: Hybrid functional, SQ, SCAN are not supported in this version of Cyclix implementation.\n" RESET);
             exit(EXIT_FAILURE);
@@ -2379,7 +2369,7 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {
             exit(EXIT_FAILURE);
         }
 
-        if (pSPARC->SOC_Flag || pSPARC->usefock || pSPARC->mGGAflag) {
+        if (pSPARC->SOC_Flag || pSPARC->usefock || pSPARC->ixc[2]) {
             if (!rank) 
                 printf(RED "ERROR: Hybrid functional, spin-orbit coupling, and SCAN are not supported in this version of SQ implementation." RESET);
             exit(EXIT_FAILURE);
@@ -3177,7 +3167,7 @@ void write_output_init(SPARC_OBJ *pSPARC) {
     }
 
     fprintf(output_fp,"***************************************************************************\n");
-    fprintf(output_fp,"*                       SPARC (version Jul 10, 2023)                      *\n");
+    fprintf(output_fp,"*                       SPARC (version Aug 01, 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);
@@ -3980,6 +3970,8 @@ Exchange Correlation
 
    Meta-GGA:     "nom"    none                           imeta=0
                  "scan"   SCAN-Meta-GGA                  imeta=1
+                 "rscan"   rSCAN-Meta-GGA                imeta=1
+                 "r2scan"  r2SCAN-Meta-GGA               imeta=1
 
    van der Waals "nov"    none                           ivdw=0
                  "vdw1"   vdW-DF1                        ivdw=1
@@ -4052,6 +4044,11 @@ void xc_decomposition(SPARC_OBJ *pSPARC)
         pSPARC->ixc[0] = 4; pSPARC->ixc[1] = 4; 
         pSPARC->ixc[2] = 1; pSPARC->ixc[3] = 0;
         pSPARC->isgradient = 1;
+    } else if (strcmpi(pSPARC->XC, "RSCAN") == 0) {
+        xc = -493494;
+        pSPARC->ixc[0] = 5; pSPARC->ixc[1] = 5; 
+        pSPARC->ixc[2] = 1; pSPARC->ixc[3] = 0;
+        pSPARC->isgradient = 1;
     } else if (strcmpi(pSPARC->XC, "R2SCAN") == 0) {
         xc = -497498;
         pSPARC->ixc[0] = 6; pSPARC->ixc[1] = 6; 

src/makefile

@@ -89,8 +89,9 @@ OBJSC = main.o initialization.o readfiles.o atomdata.o parallelization.o relax.o
         xc/vdW/vdWDF/vdWDFinitialization.o xc/vdW/vdWDF/vdWDFfinalization.o                        \
         xc/vdW/vdWDF/vdWDFnonlinearCorre.o xc/vdW/vdWDF/vdWDFreadKernel.o                 \
         xc/vdW/vdWDF/vdWDFstress.o xc/vdW/vdWDF/vdWDFparallelization.o                        \
-        xc/mgga/mGGAscan.o xc/mgga/mGGAr2scan.o xc/mgga/mGGAhamiltonianTerm.o xc/mgga/mGGAstress.o \
+        xc/mgga/mGGAscan.o xc/mgga/mGGArscan.o xc/mgga/mGGAr2scan.o xc/mgga/mGGAhamiltonianTerm.o  \
         xc/mgga/mGGAinitialization.o xc/mgga/mGGAfinalization.o xc/mgga/mGGAtauTransferTauVxc.o  \
+        xc/mgga/mGGAstress.o \
         xc/exx/exactExchange.o xc/exx/exactExchangeKpt.o xc/exx/exactExchangeInitialization.o      \
         xc/exx/exactExchangeFinalization.o xc/exx/exactExchangeStress.o                            \
         xc/exx/exactExchangePressure.o xc/exx/exactExchangeEnergyDensity.o \

src/md.c

@@ -1900,9 +1900,21 @@ void Print_fullMD(SPARC_OBJ *pSPARC, FILE *output_md, double *avgvel, double *ma
 			fprintf(output_md," %18.10E\n",maxvel[atm]);
 	#endif
 		fprintf(output_md,":MIND:\n");
-		int Nintr = pSPARC->Ntypes * (pSPARC->Ntypes + 1) / 2;
-		for(atm = 0; atm < Nintr; atm++)
-			fprintf(output_md," %18.10E\n",mindis[atm]);
+		char elemType1[8], elemType2[8]; 
+		int typeIndex, typeIndex2, pairIndex;
+		for (typeIndex = 0; typeIndex < pSPARC->Ntypes; typeIndex++) {
+			find_element(elemType1, &pSPARC->atomType[L_ATMTYPE*typeIndex]);
+			fprintf(output_md,"%s - %s: %18.10E\n", elemType1, elemType1, mindis[typeIndex]);
+		}
+		pairIndex = 0;
+		for(typeIndex = 0; typeIndex < pSPARC->Ntypes - 1; typeIndex++) {
+			find_element(elemType1, &pSPARC->atomType[L_ATMTYPE*typeIndex]);
+			for (typeIndex2 = typeIndex + 1; typeIndex2 < pSPARC->Ntypes; typeIndex2++) {
+				find_element(elemType2, &pSPARC->atomType[L_ATMTYPE*typeIndex2]);
+				fprintf(output_md,"%s - %s: %18.10E\n", elemType1, elemType2, mindis[pairIndex + pSPARC->Ntypes]);
+				pairIndex++;
+			}
+		}
 	}
 }
 
@@ -1978,7 +1990,7 @@ void MD_QOI(SPARC_OBJ *pSPARC, double *avgvel, double *maxvel, double *mindis) {
 
 	// Average and minimum distance
 	//*avgdis  = pow((pSPARC->range_x * pSPARC->range_y * pSPARC->range_z)/pSPARC->n_atom,1/3.0);
-	int atm2, ityp2, cc2;
+	int atm2, ityp2, cc2, pairIndex;
 	cc = 0;
 	for(ityp = 0; ityp < pSPARC->Ntypes; ityp++){
 	    mindis[ityp] = 1000000000.0;
@@ -1992,17 +2004,20 @@ void MD_QOI(SPARC_OBJ *pSPARC, double *avgvel, double *maxvel, double *mindis) {
 		cc += pSPARC->nAtomv[ityp];
 	}
 	cc = 0;
+	pairIndex = pSPARC->Ntypes;
 	for(ityp = 0; ityp < pSPARC->Ntypes - 1; ityp++){
 		cc2 = pSPARC->nAtomv[ityp] + cc;
 		for(ityp2 = ityp + 1; ityp2 < pSPARC->Ntypes; ityp2++){
-			mindis[pSPARC->Ntypes - 1 + ityp * (pSPARC->Ntypes - 1) + ityp2 - ityp] = 1000000000.0;
+			// mindis[pSPARC->Ntypes - 1 + ityp*(pSPARC->Ntypes-1 + pSPARC->Ntypes-1-(ityp - 1))/2 + ityp2 - ityp] = 1000000000.0;
+			mindis[pairIndex] = 1000000000.0;
 			for(atm = 0; atm < pSPARC->nAtomv[ityp]; atm++){
 				for(atm2 = 0; atm2 < pSPARC->nAtomv[ityp2]; atm2++){
 					temp = fabs(sqrt(pow(pSPARC->atom_pos[3 * (atm + cc)] - pSPARC->atom_pos[3 * (atm2 + cc2)],2.0) + pow(pSPARC->atom_pos[3 * (atm + cc) + 1] - pSPARC->atom_pos[3 * (atm2 + cc2) + 1],2.0) + pow(pSPARC->atom_pos[3 * (atm + cc) + 2] - pSPARC->atom_pos[3 * (atm2 + cc2) + 2],2.0) ));
-					if(temp < mindis[pSPARC->Ntypes - 1 + ityp * (pSPARC->Ntypes - 1) + ityp2 - ityp])
-						mindis[pSPARC->Ntypes - 1 + ityp * (pSPARC->Ntypes - 1) + ityp2 - ityp] = temp;
+					if(temp < mindis[pairIndex])
+						mindis[pairIndex] = temp;
 				}
 			}
+			pairIndex++;
 			cc2 += pSPARC->nAtomv[ityp2];
 		}
 		cc += pSPARC->nAtomv[ityp];

src/pressure.c

@@ -61,7 +61,7 @@ void Calculate_electronic_pressure(SPARC_OBJ *pSPARC) {
 
     // find exchange-correlation component of pressure
     Calculate_XC_pressure(pSPARC);
-    if ((pSPARC->mGGAflag == 1) && (pSPARC->countPotentialCalculate > 1)) { // metaGGA pressure is related to wavefunction psi directly; it needs to be computed outside of function Calculate_XC_pressure
+    if (pSPARC->ixc[2] && (pSPARC->countPotentialCalculate > 1)) { // metaGGA pressure is related to wavefunction psi directly; it needs to be computed outside of function Calculate_XC_pressure
         if (pSPARC->isGammaPoint) {
             Calculate_XC_stress_mGGA_psi_term(pSPARC); // the function is in file mgga/mgga.c
         }
@@ -204,7 +204,7 @@ void Calculate_XC_pressure(SPARC_OBJ *pSPARC) {
         // if (!rank) printf("XC pressure after d3 %9.6E\n", pSPARC->pres_xc);
     }
 
-    if (pSPARC->vdWDFFlag != 0) { // either vdW_DF1 or vdW_DF2, compute the contribution of nonlinear correlation of vdWDF on stress/pressure
+    if (pSPARC->ixc[3] != 0) { // either vdW_DF1 or vdW_DF2, compute the contribution of nonlinear correlation of vdWDF on stress/pressure
         Calculate_XC_stress_vdWDF(pSPARC); // the function is in file vdW/vdWDF/vdWDF.c
     }