Add target pressure for cell relaxation

Author: sayanbhowmik

ChangeLog

@@ -2,6 +2,14 @@
 -Date
 -Name
 -changes
+--------------
+July 17, 2025
+Name: Sayan Bhowmik
+Changes: src/include/isddft.h, src/initialization.c, src/readfiles.c, src/relax.c,
+         tests/Al_pr_relax_comp, tests/Al_pr_relax_tens, tests/SPARC_testing_script.py, 
+         tests/README.md, doc/.LaTeX/Introduction.tex, doc/.LaTeX/Optimization.tex,doc/Manual.pdf
+1. Added target pressure for cell relaxation.
+
 --------------
 June 17, 2025
 Name: Tian Tian

doc/.LaTeX/Introduction.tex

@@ -53,7 +53,7 @@
         \item \textbf{Tian Tian}: Socket communication layer, code testing, Python API
         \item \textbf{Lucas R Timmerman\footnotemark[2]}: Socket communication layer, code testing, Python API
     \end{itemize}
- 	\end{itemize}
+ \end{itemize}
 \footnotetext[1]{co-advised with Andrew J. Medford}
 \footnotetext[2]{co-advised with Phanish Suryanarayana}
 \end{frame}
@@ -378,7 +378,7 @@
   \hyperlink{RELAX_METHOD}{\texttt{RELAX\_METHOD}} $\vert$
   \hyperlink{RELAX_NITER}{\texttt{RELAX\_NITER}} $\vert$
   \hyperlink{TOL_RELAX}{\texttt{TOL\_RELAX}} $\vert$
-  \hyperlink{TOL_RELAX_CELL}{\texttt{TOL\_RELAX\_CELL}} $\vert$
+  \hyperlink{TOL_RELAX_CELL}{\texttt{TOL\_RELAX\_CELL}} $\vert$ \hyperlink{RELAX_PRESSURE}{\texttt{RELAX\_PRESSURE}} $\vert$
   \hyperlink{RELAX_MAXDILAT}{\texttt{RELAX\_MAXDILAT}} $\vert$
   \hyperlink{NLCG_SIGMA}{\texttt{NLCG\_SIGMA}} $\vert$
   \hyperlink{L_HISTORY}{\texttt{L\_HISTORY}} $\vert$

doc/.LaTeX/Optimization.tex

@@ -185,6 +185,38 @@
 \end{frame}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\begin{frame}[allowframebreaks]{\texttt{RELAX\_PRESSURE}} \label{RELAX_PRESSURE}
+\vspace*{-12pt}
+\begin{columns}
+\column{0.4\linewidth}
+\begin{block}{Type}
+Double
+\end{block}
+
+\begin{block}{Default}
+0.0
+\end{block}
+
+\column{0.4\linewidth}
+\begin{block}{Unit}
+GPa
+\end{block}
+
+\begin{block}{Example}
+\texttt{RELAX\_PRESSURE}: 4.53
+\end{block}
+\end{columns}
+
+\begin{block}{Description}
+Specifies the target outer pressure for cell volume relaxation - \hyperlink{RELAX_FLAG}{\texttt{RELAX\_FLAG}}: 2 or \hyperlink{RELAX_FLAG}{\texttt{RELAX\_FLAG}}: 3.
+\end{block}
+
+%\begin{block}{Remark}
+%\end{block}
+
+\end{frame}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \begin{frame}[allowframebreaks]{\texttt{RELAX\_MAXDILAT}} \label{RELAX_MAXDILAT}

doc/Manual.pdf

@@ -857,6 +857,7 @@ typedef struct _SPARC_OBJ{
     double qmass;          // mass parameter
     double amu2au;         // conversion factor for atomic mass unit -> atomic unit of mass
     double fs2atu;         // conversion factor for femto second -> atomic unit of time (Jiffy) 
+    double relaxPrTarget;  // Target pressure for cell relaxation in GPa
     // NPT
     int NPTscaleVecs[3];    // which lattice vector can be rescaled?
     int NPTconstraintFlag; // confinement on side length of cell. none: no length confinement (default); 1: a:b keeps unchanged; 2: a:c keeps unchanged; 
@@ -1349,6 +1350,7 @@ typedef struct _SPARC_INPUT_OBJ{
     int L_history;
     int L_autoscale;
     int L_lineopt;
+    double relaxPrTarget;  // Target pressure for cell relaxation in GPa
 
     /* Stress options*/
     int Calc_stress;

src/include/isddft.h

@@ -55,8 +55,7 @@
 #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 210
 
 
 /**
@@ -856,6 +855,7 @@ void set_defaults(SPARC_INPUT_OBJ *pSPARC_Input, SPARC_OBJ *pSPARC) {
     /* Default cell relaxation parameters*/
     pSPARC_Input->max_dilatation = 1.06;      // maximum lattice dilatation
     pSPARC_Input->TOL_RELAX_CELL = 1e-2;      // in GPa (all periodic)
+    pSPARC_Input->relaxPrTarget  = 0;         // Default relaxation pressure in 0 GPa
 
     /* Default band structure calculation parameters */
     pSPARC_Input->BandStructFlag = 0;
@@ -1575,6 +1575,7 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {
     pSPARC->OFDFT_lambda = pSPARC_Input->OFDFT_lambda;
     pSPARC->twist = pSPARC_Input->twist;
     pSPARC->is_hubbard = pSPARC_Input->is_hubbard; // DFT+U
+    pSPARC->relaxPrTarget = pSPARC_Input->relaxPrTarget; // Target pressure for cell relaxation
 
     // char type values
     strncpy(pSPARC->MDMeth , pSPARC_Input->MDMeth,sizeof(pSPARC->MDMeth));
@@ -1704,7 +1705,7 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {
     if (pSPARC->BandStructFlag == 1) {
         pSPARC->MAXIT_SCF = 1;
         pSPARC->MINIT_SCF = 1; 
-        if (pSPARC->is_hubbard) pSPARC->MAXIT_SCF = 50;
+        // if (pSPARC->is_hubbard) pSPARC->MAXIT_SCF = 50; // need to work on hubbard band
         pSPARC->PrintElecDensFlag = 0;
         pSPARC->PrintEigenFlag = 1;
     } 
@@ -3856,6 +3857,10 @@ void write_output_init(SPARC_OBJ *pSPARC) {
         }
     }
 
+    if (pSPARC->RelaxFlag > 1) {
+        fprintf(output_fp,"RELAX_PRESSURE: %.15g\n",pSPARC->relaxPrTarget);
+    }
+
     fprintf(output_fp,"CALC_STRESS: %d\n",pSPARC->Calc_stress);
     if(pSPARC->Calc_stress == 0)
         fprintf(output_fp,"CALC_PRES: %d\n",pSPARC->Calc_pres);
@@ -4355,6 +4360,7 @@ void SPARC_Input_MPI_create(MPI_Datatype *pSPARC_INPUT_MPI) {
                                          MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
                                          MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
                                          MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
+                                         MPI_DOUBLE,
                                          MPI_CHAR, MPI_CHAR, MPI_CHAR, MPI_CHAR, MPI_CHAR,
                                          MPI_CHAR, MPI_CHAR, MPI_CHAR, MPI_CHAR, MPI_CHAR};
     int blens[N_MEMBR] = {3, 3, 7,      /* int array */ 
@@ -4397,7 +4403,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,  /* double */
+                          1, 1,  /* double */
                           32, 32, 32, L_STRING, L_STRING, /* char */
                           L_STRING, L_STRING, L_STRING, L_STRING, L_STRING};
 
@@ -4611,6 +4617,8 @@ void SPARC_Input_MPI_create(MPI_Datatype *pSPARC_INPUT_MPI) {
     MPI_Get_address(&sparc_input_tmp.xi_3_SOAP, addr + i++);
     MPI_Get_address(&sparc_input_tmp.F_tol_SOAP, addr + i++);
     MPI_Get_address(&sparc_input_tmp.F_rel_scale, addr + i++);
+
+    MPI_Get_address(&sparc_input_tmp.relaxPrTarget, addr + i++);
 
 
     // char type

src/initialization.c

@@ -1014,6 +1014,9 @@ void read_input(SPARC_INPUT_OBJ *pSPARC_Input, SPARC_OBJ *pSPARC) {
             fscanf(input_fp, "%*[^\n]\n");
         } else if(strcmpi(str,"HUBBARD_FLAG:") == 0) {
             fscanf(input_fp,"%d",&pSPARC_Input->is_hubbard);
+        } else if (strcmpi(str,"RELAX_PRESSURE:") == 0) {    
+            fscanf(input_fp,"%lf",&pSPARC_Input->relaxPrTarget); // input in GPa
+            fscanf(input_fp, "%*[^\n]\n");
         } else {
             printf("\nCannot recognize input variable identifier: \"%s\"\n",str);
             exit(EXIT_FAILURE);

src/readfiles.c

@@ -1071,8 +1071,11 @@ double volrelax_constraint(SPARC_OBJ *pSPARC, double vol)
     // free memory
     free(stress_lc);
 
-    // return max pressure stress
-    return max_P_stress;
+    // return max pressure stress - target relaxation pressure
+    #ifdef DEBUG
+    if (rank == 0) printf("\n\n CELL_RELAX residual = %6.6f\n\n",-max_P_stress - pSPARC->relaxPrTarget);
+    #endif
+    return (-max_P_stress - pSPARC->relaxPrTarget);
 }
 
 

src/relax.c

@@ -0,0 +1,30 @@
+
+LATVEC_SCALE: 1 1 1  
+LATVEC:                       
+7.628635126863900 0.000000000000000 0.000000000000000 
+0.000000000000000 7.628635126863900 0.000000000000000 
+0.000000000000000 0.000000000000000 7.628635126863900
+BC: P P P                    
+
+
+EXCHANGE_CORRELATION: GGA_PBE                  
+
+
+MESH_SPACING: 0.3      
+
+KPOINT_GRID: 4 4 4          
+
+RELAX_FLAG: 2
+RELAX_PRESSURE: 1.0517108455E+01
+PRECOND_KERKER_THRESH: 0.0   
+MIXING_PARAMETER: 1.0         
+TOL_SCF: 1e-5       
+
+TOL_RELAX_CELL: 1e-4
+
+CALC_STRESS: 1
+PRINT_ATOMS: 1
+PRINT_VELS: 0
+PRINT_FORCES: 1
+PRINT_EIGEN: 0
+PRINT_DENSITY: 0

tests/Al_pr_relax_comp/high_accuracy/Al_pr_relax_comp.inpt

@@ -0,0 +1,22 @@
+#=========================
+# format of ion file
+#=========================
+# ATOM_TYPE: <atom type name>
+# N_TYPE_ATOM: <num of atoms of this type>
+# COORD:
+# <xcoord> <ycoord> <zcoord>
+# ...
+# RELAX:
+# <xrelax> <yrelax> <zrelax>
+# ...
+
+ATOM_TYPE: Al                 # atom type
+ATOMIC_MASS: 26.981539           # atomic mass (amu)
+PSEUDO_POT: ../../../psps/13_Al_3_1.9_1.9_pbe_n_v1.0.psp8  # pseudopotential file
+N_TYPE_ATOM: 4                # number of atoms of this type
+
+COORD_FRAC:                   # fractional coordinates (in lattice vector basis)
+                   0                   0                   0
+   0.500000000000000                   0   0.500000000000000
+                   0   0.500000000000000   0.500000000000000
+   0.500000000000000   0.500000000000000                   0