Fix i/o for DFT+U

ChangeLog

@@ -3,6 +3,19 @@
 -Name
 -changes
 
+--------------
+June 03, 2025
+Name: Sayan Bhowmik
+Changes: README.md, doc/.LaTeX/Introduction.tex, doc/.LaTeX/System.tex, doc/Manual.pdf,
+         src/atom/initializationAtom.c, src/atom/sparcAtom.c, src/include/isddft.h,
+         src/initialization.c, src/readfiles.c, src/xc/hubbard/hubbardForce.c,
+         src/xc/hubbard/hubbardStress.c, src/xc/locOrbRoutines.c, 
+         src/xc/hubbard/occupationMatrix.c, tests/MoO3_hubbard/, 
+         tests/NiO_spin_kpts_non_ortho_hubbard/, tests/TiCrO4_kpt_hubbard/
+1. Changed i/o format for DFT+U.
+2. Fixed minor bug in atom code.
+3. Updated doc.
+
 --------------
 June 02, 2025
 Name: Lucas Timmerman

README.md

@@ -12,6 +12,7 @@ SPARC is an open-source software package for the accurate, effcient, and scalabl
 * Spin-orbit coupling (SOC).
 * Noncollinear spin.
 * Dispersion interactions through DFT-D3, vdW-DF1, and vdW-DF2.
+* DFT+U through Dudarev's formulation.
 * Symmetry-adaption for cyclic and/or helical symmetries (Cyclix-DFT).
 * O(N) Spectral Quadrature (SQ) method.
 * On-the-fly machine-learned force field (MLFF) molecular dynamics (MD) simulations. 

doc/.LaTeX/Introduction.tex

@@ -33,7 +33,7 @@
  	    \item \textbf{Boqin Zhang}: vdW-DF, DFT-D3, meta-GGA (SCAN) \\  
  	    \item \textbf{Shashikant Kumar}: Testing framework, NLCC, MLFF \\  
             \item \textbf{Mostafa Faghih Shojaei}: SPMS table of pseudopotentials \\ 
-            \item \textbf{Sayan Bhowmik\footnotemark{}}: Atom code, DFT+U (Dudarev) \\
+            \item \textbf{Sayan Bhowmik\footnotemark}: Atom code, DFT+U (Dudarev) \\
             \item \textbf{Swarnava Ghosh}: Preliminary development \\
             \item \textbf{Deepa Phanish}: Initial development
  	\end{itemize}
@@ -282,6 +282,7 @@
   \hyperlink{COORD_FRAC}{\texttt{COORD\_FRAC}} $\vert$
   \hyperlink{RELAX}{\texttt{RELAX}} $\vert$
   \hyperlink{SPIN}{\texttt{SPIN}} $\vert$
+  \hyperlink{HUBBARD_FLAG}{\texttt{HUBBARD\_FLAG}} $\vert$
   \hyperlink{HUBBARD}{\texttt{HUBBARD}} $\vert$
   \hyperlink{U_ATOM_TYPE}{\texttt{U\_ATOM\_TYPE}} $\vert$
   \hyperlink{U_VAL}{\texttt{U\_VAL}} 

doc/.LaTeX/System.tex

@@ -777,6 +777,40 @@
 \end{frame}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\begin{frame}[allowframebreaks]{\texttt{HUBBARD\_FLAG}} \label{HUBBARD_FLAG}
+\vspace*{-12pt}
+\begin{columns}
+\column{0.4\linewidth}
+\begin{block}{Type}
+0 or 1
+\end{block}
+
+\begin{block}{Default}
+\texttt{0}
+\end{block}
+
+\column{0.4\linewidth}
+\begin{block}{Unit}
+No unit
+\end{block}
+
+\begin{block}{Example}
+\texttt{HUBBARD\_FLAG}: \texttt{1}
+\end{block}
+\end{columns}
+
+\begin{block}{Description}
+Flag for adding hubbard correction on top of the \hyperlink{EXCHANGE_CORRELATION}{\texttt{EXCHANGE\_CORRELATION}}.
+\end{block}
+
+\begin{block}{Remark}
+Usually used along with LDA, GGA or SCAN type functionals.
+\end{block}
+
+\end{frame}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \begin{frame}[allowframebreaks,c]{} \label{System:ion}
@@ -1053,13 +1087,17 @@
 
 \begin{block}{Example}
 \texttt{HUBBARD: \\}
-\texttt{U\_ATOM\_TYPE: Ni \\}
-\texttt{U\_VAL:} 0 0 0.05 0 \\
+% \texttt{U\_ATOM\_TYPE: Ni \\}
+% \texttt{U\_VAL:} 0 0 0.05 0 \\
 \end{block}
 \end{columns}
 
 \begin{block}{Description}
-Triggers a DFT+U calculation on top of the \hyperlink{EXCHANGE_CORRELATION}{\texttt{EXCHANGE\_CORRELATION}} specified in the \texttt{.inpt} file. Must be followed by specifying the atoms on which a U correction is desired along with value of U as per Dudarev's scheme.
+Start of the DFT+U parameter block. Must be followed by specifying the atoms on which a U correction is desired along with value of U as per Dudarev's scheme.
+\end{block}
+
+\begin{block}{Remark}
+\hyperlink{HUBBARD_FLAG}{HUBBARD\_FLAG} must be specified to 1 in \texttt{.inpt} file.
 \end{block}
 
 \end{frame}

src/atom/initializationAtom.c

@@ -249,7 +249,7 @@ void copy_PSP_atom(SPARC_ATOM_OBJ *pSPARC_ATOM, SPARC_OBJ *pSPARC, int ityp) {
     for (int l = 0; l <= pSPARC_ATOM->psd->lmax; l++) {
         pSPARC_ATOM->psd->rc[l] = pSPARC->psd[ityp].rc[l];
         pSPARC_ATOM->psd->ppl[l] = pSPARC->psd[ityp].ppl[l];
-        lpos[l+1] = lpos[l] + pSPARC->psd->ppl[l];
+        lpos[l+1] = lpos[l] + pSPARC_ATOM->psd->ppl[l];
     }
     int nproj = lpos[lmax+1];
 

src/atom/sparcAtom.c

@@ -34,12 +34,12 @@ void sparc_atom(SPARC_OBJ *pSPARC, int ityp, SPARC_INPUT_OBJ *pSPARC_Input) {
     electronicGroundState_atom(&pSPARC_ATOM);
 
     // Printing
-    if (!pSPARC->is_hubbard){
+    if (!pSPARC_Input->is_hubbard){
         printResultsAtom(&pSPARC_ATOM);
     } 
 
     // Copy solution
-    if (pSPARC->is_hubbard) copyAtomSolution(&pSPARC_ATOM, pSPARC, ityp);
+    if (pSPARC_Input->is_hubbard) copyAtomSolution(&pSPARC_ATOM, pSPARC, ityp);
 
     // Free memory
     Finalize_Atom(&pSPARC_ATOM);

src/include/isddft.h

@@ -1471,6 +1471,9 @@ typedef struct _SPARC_INPUT_OBJ{
     double exx_frac;        // hybrid mixing coefficient    
     int ExxMethod;          // method for solving poissons equation, kronecker product way or FFT
 
+    /* DFT+U */
+    int is_hubbard;         // DFT+U flag
+
     /* SQ methods */
     int sqAmbientFlag;             // Flag of SQ method
     // int SQ_highT_gauss_mem;       // Memory option for gauss quadrature 

src/initialization.c

@@ -55,8 +55,8 @@
 #define min(x,y) ((x)<(y)?(x):(y))
 #define max(x,y) ((x)>(y)?(x):(y))
 
-//#define N_MEMBR 209
-#define N_MEMBR 208
+#define N_MEMBR 209
+// #define N_MEMBR 208
 
 
 /**
@@ -300,8 +300,8 @@ void Initialize(SPARC_OBJ *pSPARC, int argc, char *argv[]) {
         t1 = MPI_Wtime();
 #endif
 
-        // broadcast hubbard flag
-        MPI_Bcast(&pSPARC->is_hubbard, 1, MPI_INT, 0, MPI_COMM_WORLD);
+        // // broadcast hubbard flag
+        // MPI_Bcast(&pSPARC->is_hubbard, 1, MPI_INT, 0, MPI_COMM_WORLD);
 
         // broadcast Ntypes read from ion file
         MPI_Ibcast(&pSPARC->Ntypes, 1, MPI_INT, 0, MPI_COMM_WORLD, &req);
@@ -341,8 +341,8 @@ void Initialize(SPARC_OBJ *pSPARC, int argc, char *argv[]) {
         t2 = MPI_Wtime();
         if (rank == 0) printf("Broadcasting the input parameters took %.3f ms\n",(t2-t1)*1000);
 #endif
-        // broadcast hubbard flag
-        MPI_Bcast(&pSPARC->is_hubbard, 1, MPI_INT, 0, MPI_COMM_WORLD);
+        // // broadcast hubbard flag
+        // MPI_Bcast(&pSPARC->is_hubbard, 1, MPI_INT, 0, MPI_COMM_WORLD);
 
         // broadcast Ntypes read from ion file
         MPI_Ibcast(&pSPARC->Ntypes, 1, MPI_INT, 0, MPI_COMM_WORLD, &req);
@@ -967,6 +967,9 @@ void set_defaults(SPARC_INPUT_OBJ *pSPARC_Input, SPARC_OBJ *pSPARC) {
     // read in initial density
     pSPARC_Input->readInitDens = 0;
 
+    // DFT+U
+    pSPARC_Input->is_hubbard = 0;
+
     /* Default socket options
        Note to future developers: please keep the USE_SOCKET macro
        as the LAST PART of the initialization function!!
@@ -1571,6 +1574,7 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {
     pSPARC->OFDFT_tol = pSPARC_Input->OFDFT_tol;
     pSPARC->OFDFT_lambda = pSPARC_Input->OFDFT_lambda;
     pSPARC->twist = pSPARC_Input->twist;
+    pSPARC->is_hubbard = pSPARC_Input->is_hubbard; // DFT+U
 
     // char type values
     strncpy(pSPARC->MDMeth , pSPARC_Input->MDMeth,sizeof(pSPARC->MDMeth));
@@ -3777,6 +3781,8 @@ void write_output_init(SPARC_OBJ *pSPARC) {
         fprintf(output_fp,"EXX_RANGE_FOCK: %.6f\n", pSPARC->hyb_range_fock);
         fprintf(output_fp,"EXX_RANGE_PBE: %.6f\n", pSPARC->hyb_range_pbe);
     }
+    // DFT+U
+    fprintf(output_fp, "HUBBARD_FLAG: %d\n", pSPARC->is_hubbard);
     if (pSPARC->sqAmbientFlag == 1 || pSPARC->sqHighTFlag == 1) {
         if (pSPARC->sqAmbientFlag) fprintf(output_fp,"SQ_AMBIENT_FLAG: %d\n", pSPARC->sqAmbientFlag);
         if (pSPARC->sqHighTFlag) fprintf(output_fp,"SQ_HIGHT_FLAG: %d\n", pSPARC->sqHighTFlag);
@@ -4056,6 +4062,7 @@ void write_output_init(SPARC_OBJ *pSPARC) {
                 pSPARC->stress_rel_scale[4], pSPARC->stress_rel_scale[5]);
         fprintf(output_fp, "MLFF_DFT_FQ: %d\n", pSPARC->MLFF_DFT_fq);
     }
+
 
     fprintf(output_fp,"VERBOSITY: %d\n",pSPARC->Verbosity);
     fprintf(output_fp,"PRINT_FORCES: %d\n",pSPARC->PrintForceFlag);
@@ -4329,7 +4336,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, /* int array */
+                                         MPI_INT, MPI_INT, MPI_INT, MPI_INT, /* int array */
 
                                          MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
                                          MPI_DOUBLE, MPI_DOUBLE,
@@ -4373,7 +4380,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, /* int */ 
+                          1, 1, 1, 1, /* int */ 
                           9, 3, L_QMASS, L_kpoint, L_kpoint,
                           L_kpoint, 6, /* double array */
                           1, 1, 1, 1, 1, 
@@ -4520,6 +4527,7 @@ void SPARC_Input_MPI_create(MPI_Datatype *pSPARC_INPUT_MPI) {
     MPI_Get_address(&sparc_input_tmp.N_rgrid_MLFF, addr + i++);
     MPI_Get_address(&sparc_input_tmp.MLFF_DFT_fq, addr + i++);
     MPI_Get_address(&sparc_input_tmp.REFERENCE_CUTOFF_FAC, addr + i++);
+    MPI_Get_address(&sparc_input_tmp.is_hubbard, addr + i++);
 
     // double array type
     MPI_Get_address(&sparc_input_tmp.LatVec, addr + i++);

src/readfiles.c

@@ -1012,6 +1012,8 @@ void read_input(SPARC_INPUT_OBJ *pSPARC_Input, SPARC_OBJ *pSPARC) {
         } else if(strcmpi(str,"PRINT_ENERGY_DENSITY:") == 0) {
             fscanf(input_fp,"%d",&pSPARC_Input->PrintEnergyDensFlag);
             fscanf(input_fp, "%*[^\n]\n");
+        } else if(strcmpi(str,"HUBBARD_FLAG:") == 0) {
+            fscanf(input_fp,"%d",&pSPARC_Input->is_hubbard);
         } else {
             printf("\nCannot recognize input variable identifier: \"%s\"\n",str);
             exit(EXIT_FAILURE);
@@ -1848,7 +1850,7 @@ void read_ion(SPARC_INPUT_OBJ *pSPARC_Input, SPARC_OBJ *pSPARC) {
 
     //////////////////////////// HUBBARD BLOCK ///////////////////////////////////
     // @author: Sayan Bhowmik (sbhowmik9@gatech.edu)
-    pSPARC->is_hubbard = 0;
+    // pSPARC->is_hubbard = 0;
     fseek(ion_fp, 0L, SEEK_SET);
 
     /* Identify atom types for which radial solve is desired */
@@ -1869,9 +1871,13 @@ void read_ion(SPARC_INPUT_OBJ *pSPARC_Input, SPARC_OBJ *pSPARC) {
 #ifdef DEBUG
             printf("Found Hubbard block.\n");
 #endif
+            if (pSPARC_Input->is_hubbard == 0) {
+                printf("Please specify HUBBARD_FALG: 1 in the .inpt file.\n");
+                exit(EXIT_FAILURE);
+            }
             pSPARC->atomSolvSpinFlag = 0;
             atm_flag = 1;
-            pSPARC->is_hubbard = 1;
+            // pSPARC->is_hubbard = 1;
             pSPARC->AtmU = (HUB_DUDAREV_OBJ *)malloc(pSPARC->Ntypes*sizeof(HUB_DUDAREV_OBJ)); // initialise the AtmU object
             fscanf(ion_fp, "%*[^\n]\n"); 
         } else if (strcmpi(str, "U_ATOM_TYPE:") == 0) {
@@ -1896,7 +1902,8 @@ void read_ion(SPARC_INPUT_OBJ *pSPARC_Input, SPARC_OBJ *pSPARC) {
             } else {
                 typcnt++;
                 if (pSPARC->atom_solve_flag != NULL) {
-                    pSPARC->atom_solve_flag[typcnt] = 1;
+                    // pSPARC->atom_solve_flag[typcnt] = 1;
+                    pSPARC->atom_solve_flag[id_OG] = 1;
                 } else {
                     printf("Check if you have indicated start of the HUBBARD block...\n");
                     printf("Include the following line to indicate start of the block\n");
@@ -1907,8 +1914,10 @@ void read_ion(SPARC_INPUT_OBJ *pSPARC_Input, SPARC_OBJ *pSPARC) {
 
             free(local_str);
         } else if(strcmpi(str, "U_VAL:") == 0) {
-            Uval_count = fscanf(ion_fp,"%lf %lf %lf %lf",&pSPARC->AtmU[typcnt].U[0], &pSPARC->AtmU[typcnt].U[1],
-                &pSPARC->AtmU[typcnt].U[2], &pSPARC->AtmU[typcnt].U[3]);
+            // Uval_count = fscanf(ion_fp,"%lf %lf %lf %lf",&pSPARC->AtmU[typcnt].U[0], &pSPARC->AtmU[typcnt].U[1],
+            //     &pSPARC->AtmU[typcnt].U[2], &pSPARC->AtmU[typcnt].U[3]);
+            Uval_count = fscanf(ion_fp,"%lf %lf %lf %lf",&pSPARC->AtmU[id_OG].U[0], &pSPARC->AtmU[id_OG].U[1],
+                &pSPARC->AtmU[id_OG].U[2], &pSPARC->AtmU[id_OG].U[3]);
             if (Uval_count != 4) {
                 printf("\nError: U_VAL must have exactly 4 values corresponding to s p d f orbitals. Found %d values.\n", Uval_count);
                 exit(EXIT_FAILURE);
@@ -1923,7 +1932,7 @@ void read_ion(SPARC_INPUT_OBJ *pSPARC_Input, SPARC_OBJ *pSPARC) {
     }
 
     // checks
-    if ((pSPARC->is_hubbard == 1) && (typcnt+1 == 0)) {
+    if ((pSPARC_Input->is_hubbard == 1) && (typcnt+1 == 0)) {
         printf("\nPlease specify atoms where you want to apply hubbard corrections.\n");
         exit(EXIT_FAILURE);
     }

src/xc/hubbard/hubbardForce.c

@@ -227,7 +227,7 @@ void Compute_Integral_Orb_Dpsi(SPARC_OBJ *pSPARC, double *dpsi, double *beta) {
  * @brief Calculate forces in gamma point using integrals.
  */
 void Compute_force_hubbard_by_integrals(SPARC_OBJ *pSPARC, double *force_hub, double *alpha) {
-    int n, ityp, iat, ncol, atom_index, atom_index2, count;
+    int n, ityp, iat, ncol, atom_index, atom_index2, first_id, count;
     int spinor, Nspinor, hub_sp;
     double *pre_fac;
     int angnum;
@@ -283,10 +283,15 @@ void Compute_force_hubbard_by_integrals(SPARC_OBJ *pSPARC, double *force_hub, do
     for (spinor = 0; spinor < Nspinor; spinor++) {
         int spinorshift = pSPARC->IP_displ_U[atm_idx] * ncol * spinor;
 
-        atom_index = 0; 
+        atom_index = 0; first_id = 0;
 
         for (ityp = 0; ityp < pSPARC->Ntypes; ityp++) {
-            if (!pSPARC->atom_solve_flag[ityp]) continue;
+            if (!pSPARC->atom_solve_flag[ityp]) {
+                if (first_id == 0) atom_index += pSPARC->nAtomv[ityp];
+                continue;
+            }
+
+            if (first_id == 0) first_id++;
 
             angnum = pSPARC->AtmU[ityp].angnum;
             for (iat = 0; iat < pSPARC->nAtomv[ityp]; iat++) {
@@ -315,14 +320,17 @@ void Compute_force_hubbard_by_integrals(SPARC_OBJ *pSPARC, double *force_hub, do
             hub_sp = spinor;
         }
 
-        atom_index = 0; atom_index2 = 0;
+        atom_index = 0; atom_index2 = 0, first_id = 0;
 
         for (ityp = 0; ityp < pSPARC->Ntypes; ityp++) {
             if (!pSPARC->atom_solve_flag[ityp]) {
-                continue;
                 atom_index += pSPARC->nAtomv[ityp];
+                if (first_id == 0) atom_index2 = atom_index;
+                continue;
             }
 
+            if (first_id == 0) first_id++;
+
             angnum = pSPARC->AtmU[ityp].angnum;
 
             for (iat = 0; iat < pSPARC->nAtomv[ityp]; iat++) {
@@ -663,7 +671,7 @@ void Compute_Integral_Orb_Dpsi_kpt(SPARC_OBJ *pSPARC, double _Complex *dpsi, dou
  * @brief Calculate k-point DFT+U forces using integral.
  */
 void Compute_force_hubbard_by_integrals_kpt(SPARC_OBJ *pSPARC, double *force_hub, double _Complex *alpha) {
-    int k, n, ityp, iat, ncol, atom_index, atom_index2, count, Nk;
+    int k, n, ityp, iat, ncol, atom_index, atom_index2, first_id, count, Nk;
     int spinor, Nspinor, hub_sp, *IP_displ_U;
     ncol = pSPARC->Nband_bandcomm; // number of bands assigned
     Nk = pSPARC->Nkpts_kptcomm;
@@ -721,10 +729,15 @@ void Compute_force_hubbard_by_integrals_kpt(SPARC_OBJ *pSPARC, double *force_hub
         for (spinor = 0; spinor < Nspinor; spinor++) {
             int spinorshift = pSPARC->IP_displ_U[atm_idx] * ncol * Nspinor * k + pSPARC->IP_displ_U[atm_idx] * ncol * spinor;
 
-            atom_index = 0;
+            atom_index = 0; first_id = 0;
 
             for (ityp = 0; ityp < pSPARC->Ntypes; ityp++) {
-                if (!pSPARC->atom_solve_flag[ityp]) continue;
+                if (!pSPARC->atom_solve_flag[ityp]) {
+                    if (first_id == 0) atom_index += pSPARC->nAtomv[ityp];
+                    continue;
+                }
+
+                if (first_id == 0) first_id++;
 
                 angnum = pSPARC->AtmU[ityp].angnum;
                 for (iat = 0; iat < pSPARC->nAtomv[ityp]; iat++) {
@@ -756,14 +769,17 @@ void Compute_force_hubbard_by_integrals_kpt(SPARC_OBJ *pSPARC, double *force_hub
                 hub_sp = spinor;
             }
 
-            atom_index = 0; atom_index2 = 0;
+            atom_index = 0; atom_index2 = 0, first_id = 0;
 
             for (ityp = 0; ityp < pSPARC->Ntypes; ityp++) {
                 if (!pSPARC->atom_solve_flag[ityp]) {
-                    continue;
                     atom_index += pSPARC->nAtomv[ityp];
+                    if (first_id == 0) atom_index2 = atom_index;
+                    continue;
                 }
 
+                if (first_id == 0) first_id++;
+
                 angnum = pSPARC->AtmU[ityp].angnum;
 
                 for (iat = 0; iat < pSPARC->nAtomv[ityp]; iat++) {