Skip to content

Cyclix Stress bug fixed #252

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
abhirajsharma wants to merge 1 commit into from
Closed

Cyclix Stress bug fixed #252

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions ChangeLog
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,14 @@
-Name
-changes
--------------
September 25, 2025
Name: Rajat kumar, Abhiraj Sharma
Changes: src/cyclix/cyclix_stress.c, src/cyclix/cyclix_stress.h, src/stress.c, tests/cyclix
1. Bug fixes in stress in cyclix due to more than 1 kpts/processor as well as spin calculation
2. Bug fix in stress units in .out file as well as Ha/Bohr to GPa conversion
3. Fixed reference results both standard as well as high accuracy in cyclix


September 04, 2025
Name: Sayan Bhowmik
Changes: (src/md.c)
Expand Down
44 changes: 20 additions & 24 deletions src/cyclix/cyclix_stress.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -853,17 +853,17 @@ void Compute_stress_tensor_kinetic_cyclix(SPARC_OBJ *pSPARC, double *dpsi, doubl

double temp_k = 0.0;
for(n = 0; n < ncol; n++){
double dpsi3_dpsi3 = 0.0;
for (spinor = 0; spinor < Nspinor; spinor++) {
double dpsi3_dpsi3 = 0.0;
double *dpsi_ptr = dpsi + n * DMndsp + spinor * DMnd; // dpsi_1
for(i = 0; i < DMnd; i++){
dpsi3_dpsi3 += *(dpsi_ptr + i) * *(dpsi_ptr + i) * pSPARC->Intgwt_psi[i];
}
}
double *occ = pSPARC->occ;
if (pSPARC->spin_typ == 1) occ += spinor * Ns;
double g_nk = occ[n + pSPARC->band_start_indx];
temp_k += dpsi3_dpsi3 * g_nk;
double *occ = pSPARC->occ;
if (pSPARC->spin_typ == 1) occ += spinor * Ns;
double g_nk = occ[n + pSPARC->band_start_indx];
temp_k += dpsi3_dpsi3 * g_nk;
}
}
*stress_k = - pSPARC->occfac * temp_k;
}
Expand All @@ -882,7 +882,7 @@ void Compute_stress_tensor_nloc_by_integrals_cyclix(SPARC_OBJ *pSPARC, double *s
IP_displ = pSPARC->IP_displ;

double *beta3_x3 = alpha + IP_displ[pSPARC->n_atom]*ncol*Nspinor;
double stress_nl_ = 0;
double stress_nl_ = 0.0;
count = 0;
for (spinor = 0; spinor < Nspinor; spinor++) {
for (ityp = 0; ityp < pSPARC->Ntypes; ityp++) {
Expand Down Expand Up @@ -984,9 +984,9 @@ void Calculate_nonlocal_kinetic_stress_kpt_cyclix(SPARC_OBJ *pSPARC)
beta3 = alpha_so2 + pSPARC->IP_displ_SOC[pSPARC->n_atom] * ncol * Nspinor * (Nk + kpt);
Compute_Integral_Chi_XmRjp_beta_Dpsi_kpt_cyclix(pSPARC, pSPARC->Yorb_kpt, beta3, kpt, "SO2");
}
}

Compute_stress_tensor_kinetic_kpt_cyclix(pSPARC, pSPARC->Yorb_kpt, &stress_k);
Compute_stress_tensor_kinetic_kpt_cyclix(pSPARC, pSPARC->Yorb_kpt, &stress_k, kpt);
}

if (pSPARC->npNd > 1) {
MPI_Allreduce(MPI_IN_PLACE, alpha, pSPARC->IP_displ[pSPARC->n_atom] * ncol * Nk * 2 * Nspinor, MPI_DOUBLE_COMPLEX, MPI_SUM, pSPARC->dmcomm);
Expand Down Expand Up @@ -1137,7 +1137,7 @@ void Compute_Integral_Chi_XmRjp_beta_Dpsi_kpt_cyclix(SPARC_OBJ *pSPARC, double _
/**
* @brief Calculate kinetic stress tensor
*/
void Compute_stress_tensor_kinetic_kpt_cyclix(SPARC_OBJ *pSPARC, double _Complex *dpsi, double *stress_k)
void Compute_stress_tensor_kinetic_kpt_cyclix(SPARC_OBJ *pSPARC, double _Complex *dpsi, double *stress_k, int kpt)
{
int ncol, DMnd, Nspinor, DMndsp, Nk, Ns;
ncol = pSPARC->Nband_bandcomm; // number of bands assigned
Expand All @@ -1147,26 +1147,22 @@ void Compute_stress_tensor_kinetic_kpt_cyclix(SPARC_OBJ *pSPARC, double _Complex
Nk = pSPARC->Nkpts_kptcomm;
Ns = pSPARC->Nstates;

int kpt, n, spinor, i;
double stress_k_ = 0;
for(kpt = 0; kpt < Nk; kpt++) {
double temp_k = 0.0;
for(n = 0; n < ncol; n++){
int n, spinor, i;
double temp_k = 0.0;
for(n = 0; n < ncol; n++){
for (spinor = 0; spinor < Nspinor; spinor++) {
double dpsi3_dpsi3 = 0.0;
for (spinor = 0; spinor < Nspinor; spinor++) {
double _Complex *dpsi_ptr = dpsi + n * DMndsp + spinor * DMnd; // dpsi_1
for(i = 0; i < DMnd; i++){
dpsi3_dpsi3 += (creal(*(dpsi_ptr + i)) * creal(*(dpsi_ptr + i)) + cimag(*(dpsi_ptr + i)) * cimag(*(dpsi_ptr + i))) * pSPARC->Intgwt_psi[i];
}
double _Complex *dpsi_ptr = dpsi + n * DMndsp + spinor * DMnd; // dpsi_1
for(i = 0; i < DMnd; i++){
dpsi3_dpsi3 += (creal(*(dpsi_ptr + i)) * creal(*(dpsi_ptr + i)) + cimag(*(dpsi_ptr + i)) * cimag(*(dpsi_ptr + i))) * pSPARC->Intgwt_psi[i];
}
double *occ = pSPARC->occ + kpt*Ns;
if (pSPARC->spin_typ == 1) occ += spinor * Nk * Ns;
double g_nk = occ[n + pSPARC->band_start_indx];
temp_k += dpsi3_dpsi3 * g_nk;
}
stress_k_ -= pSPARC->occfac * pSPARC->kptWts_loc[kpt] / pSPARC->Nkpts * temp_k;
}
*stress_k = stress_k_;
*stress_k -= pSPARC->occfac * pSPARC->kptWts_loc[kpt] / pSPARC->Nkpts * temp_k;
}


Expand All @@ -1175,7 +1171,7 @@ void Compute_stress_tensor_nloc_by_integrals_kpt_cyclix(SPARC_OBJ *pSPARC, doubl
int k, n, np, ldispl, ityp, iat, ncol, Ns;
int count, l, m, lmax, Nk, spinor, Nspinor;
int l_start, mexclude, ppl, *IP_displ;
double g_nk, kptwt, alpha_r, alpha_i, SJ, temp2_s, scaled_gamma_Jl = 0;
double g_nk, kptwt, alpha_r, alpha_i, SJ, temp2_s, scaled_gamma_Jl = 0.0;
ncol = pSPARC->Nband_bandcomm; // number of bands assigned
Ns = pSPARC->Nstates;
Nk = pSPARC->Nkpts_kptcomm;
Expand Down Expand Up @@ -1233,4 +1229,4 @@ void Compute_stress_tensor_nloc_by_integrals_kpt_cyclix(SPARC_OBJ *pSPARC, doubl
}
}
}
}
}
2 changes: 1 addition & 1 deletion src/cyclix/include/cyclix_stress.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ void Dpseudopot_cyclix_z(SPARC_OBJ *pSPARC, double *VJ, int FDn, int ishift_p, i
* @brief Calculate kinetic stress tensor
*/
void Compute_stress_tensor_kinetic_cyclix(SPARC_OBJ *pSPARC, double *dpsi, double *stress_k);
void Compute_stress_tensor_kinetic_kpt_cyclix(SPARC_OBJ *pSPARC, double _Complex *dpsi, double *stress_k);
void Compute_stress_tensor_kinetic_kpt_cyclix(SPARC_OBJ *pSPARC, double _Complex *dpsi, double *stress_k, int kpt);

/**
* @brief Calculate nonlocal stress components
Expand Down
2 changes: 1 addition & 1 deletion src/electronicGroundState.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -284,7 +284,7 @@ void Calculate_electronicGroundState(SPARC_OBJ *pSPARC) {
double cellsizes[3] = {pSPARC->range_x, pSPARC->range_y, pSPARC->range_z};
int BCs[3] = {pSPARC->BCx, pSPARC->BCy, pSPARC->BCz};
char stressUnit[16];
double maxSGPa = convertStressToGPa(maxS, cellsizes, BCs, stressUnit);
double maxSGPa = convertStressToGPa(pSPARC, maxS, cellsizes, BCs, stressUnit);
fprintf(output_fp,"Maximum stress :%18.10E (%s)\n",maxS,stressUnit);
fprintf(output_fp,"Maximum stress equiv. to periodic :%18.10E (GPa)\n",maxSGPa);
}
Expand Down
2 changes: 1 addition & 1 deletion src/include/stress.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ void PrintStress (SPARC_OBJ *pSPARC, double *stress, FILE *fp);
* @param origUnit (OUTPUT) Unit for the original stress.
* @return double Stress component in GPa.
*/
double convertStressToGPa(double Stress, double cellsizes[3], int BCs[3], char origUnit[16]);
double convertStressToGPa(SPARC_OBJ *pSPARC, double Stress, double cellsizes[3], int BCs[3], char origUnit[16]);


#endif // STRESS_H
2 changes: 1 addition & 1 deletion src/initialization.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3722,7 +3722,7 @@ void write_output_init(SPARC_OBJ *pSPARC) {
}

fprintf(output_fp,"***************************************************************************\n");
fprintf(output_fp,"* SPARC (version September 04, 2025) *\n");
fprintf(output_fp,"* SPARC (version September 25, 2025) *\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);
Expand Down
2 changes: 1 addition & 1 deletion src/makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -159,4 +159,4 @@ clean:
rm -f $(OBJSC) $(LIBBASE)
rm -f socket/*.o
test: ../tests/SPARC_testing_script.py
cd ../tests; python SPARC_testing_script.py
cd ../tests; python SPARC_testing_script.py
11 changes: 7 additions & 4 deletions src/stress.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2113,22 +2113,25 @@ void PrintStress (SPARC_OBJ *pSPARC, double *stress, FILE *fp) {
* @param origUnit (OUTPUT) Unit for the original stress.
* @return double Stress component in GPa.
*/
double convertStressToGPa(double Stress, double cellsizes[3], int BCs[3], char origUnit[16]) {
double convertStressToGPa(SPARC_OBJ *pSPARC, double Stress, double cellsizes[3], int BCs[3], char origUnit[16]) {
double scale = 1.0; // scale the stress by the dimensions of the dirichlet directions
int nperiods = 3;
if (BCs[0] == 1) { scale *= cellsizes[0]; nperiods--; }
if (BCs[1] == 1) { scale *= cellsizes[1]; nperiods--; }
if (BCs[2] == 1) { scale *= cellsizes[2]; nperiods--; }

// scale and then convert to GPa
double StressGPa = Stress / scale * CONST_HA_BOHR3_GPA;
double StressGPa;
if (pSPARC->CyclixFlag)
StressGPa = Stress*CONST_HA_BOHR3_GPA/(M_PI * ( pow(pSPARC->xout,2.0) - pow(pSPARC->xin,2.0) ) );
else
StressGPa = Stress / scale * CONST_HA_BOHR3_GPA;

// find the original unit
if (nperiods == 1)
if (nperiods == 1 || pSPARC->CyclixFlag)
snprintf(origUnit, 16, "Ha/Bohr");
else
snprintf(origUnit, 16, "Ha/Bohr**%d", nperiods);

return StressGPa;
}

Expand Down
8 changes: 4 additions & 4 deletions tests/SPARC_testing_script.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,12 +12,12 @@
import math

# Other parameters to run the test (can be changed by the user)
nprocs_tests = 24 # In default tests are run with 24 processors per node
nnodes_tests = 2 # In default tests are run with 1 node
npbs = 10 # By default (number of script files the tests are distributed to)
nprocs_tests = 48 # In default tests are run with 24 processors per node
nnodes_tests = 1 # In default tests are run with 1 node
npbs = 3 # By default (number of script files the tests are distributed to)
launch_cluster_extension = ".sbatch" # extension of the file used to launch the jobs on the cluster by default it is .sbatch
command_launch_extension = "sbatch" # Command to launch the script to ask for resources on the cluster (example: qsub launch.pbs)
MPI_command = "srun" # MPI command to run the executable on the given cluster
MPI_command = "mpirun -np 48" # MPI command to run the executable on the given cluster



Expand Down
88 changes: 45 additions & 43 deletions tests/cyclix/FeCl2_cyclix_spin/high_accuracy/FeCl2_cyclix_spin.refout
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
***************************************************************************
* SPARC (version June 24, 2024) *
* SPARC (version September 04, 2025) *
* Copyright (c) 2020 Material Physics & Mechanics Group, Georgia Tech *
* Distributed under GNU General Public License 3 (GPL) *
* Start time: Mon Jun 24 15:30:15 2024 *
* Start time: Wed Sep 24 20:42:08 2025 *
***************************************************************************
Input parameters
***************************************************************************
Expand All @@ -17,6 +17,7 @@ SPIN_TYP: 1
ELEC_TEMP_TYPE: Fermi-Dirac
SMEARING: 0.001
EXCHANGE_CORRELATION: GGA_PBE
HUBBARD_FLAG: 0
NSTATES: 41
CHEB_DEGREE: 100
CHEFSI_BOUND_FLAG: 0
Expand Down Expand Up @@ -47,7 +48,7 @@ PRINT_ATOMS: 1
PRINT_EIGEN: 0
PRINT_DENSITY: 0
PRINT_ENERGY_DENSITY: 0
OUTPUT_FILE: FeCl2_cyclix_spin/temp_run/FeCl2_cyclix_spin
OUTPUT_FILE: cyclix/FeCl2_cyclix_spin/temp_run/FeCl2_cyclix_spin
***************************************************************************
Cell
***************************************************************************
Expand All @@ -58,18 +59,18 @@ Density: 1.2141297401E-01 (amu/Bohr^3), 1.3605383749E+00 (g/cc)
***************************************************************************
NP_SPIN_PARAL: 2
NP_KPOINT_PARAL: 1
NP_BAND_PARAL: 6
NP_DOMAIN_PARAL: 4 1 2
NP_DOMAIN_PHI_PARAL: 12 2 4
NP_BAND_PARAL: 7
NP_DOMAIN_PARAL: 17 1 1
NP_DOMAIN_PHI_PARAL: 16 3 5
EIG_SERIAL_MAXNS: 1500
***************************************************************************
Initialization
***************************************************************************
Number of processors : 96
Number of processors : 240
Mesh spacing in x-direction : 0.120234 (Bohr)
Mesh spacing in y-direction : 0.00296377 (Bohr)
Mesh spacing in z-direction : 0.119703 (Bohr)
Output printed to : FeCl2_cyclix_spin/temp_run/FeCl2_cyclix_spin.out
Output printed to : cyclix/FeCl2_cyclix_spin/temp_run/FeCl2_cyclix_spin.out
Total number of atom types : 2
Total number of atoms : 6
Total number of electrons : 60
Expand All @@ -84,53 +85,54 @@ Atomic mass : 35.4515
Pseudocharge radii of atom type 2 : 6.61 0.15 6.58 (x, y, z dir)
Number of atoms of type 2 : 4
Estimated total memory usage : 5.82 GB
Estimated memory per processor : 62.09 MB
Estimated memory per processor : 24.84 MB
========================================================================================
Self Consistent Field (SCF#1)
========================================================================================
Iteration Free Energy (Ha/atom) Magnetization SCF Error Timing (sec)
1 -4.9866567407E+01 8.0000E+00 2.275E-01 18.758
2 -4.9871693480E+01 8.0000E+00 2.483E-01 6.083
3 -4.9931035419E+01 8.0000E+00 8.737E-02 5.879
4 -4.9938547782E+01 8.0000E+00 3.884E-02 5.766
5 -4.9935561172E+01 8.0000E+00 4.924E-02 5.708
6 -4.9941254720E+01 8.0000E+00 1.354E-02 5.598
7 -4.9941805651E+01 8.0000E+00 7.158E-03 5.669
8 -4.9941883346E+01 8.0000E+00 2.936E-03 5.329
9 -4.9941900215E+01 8.0000E+00 2.487E-03 5.001
10 -4.9941910449E+01 8.0000E+00 1.106E-03 5.320
11 -4.9941912265E+01 8.0000E+00 4.472E-04 8.508
12 -4.9941912597E+01 8.0000E+00 1.529E-04 5.152
13 -4.9941912656E+01 8.0000E+00 8.555E-05 4.903
14 -4.9941912686E+01 8.0000E+00 4.390E-05 4.870
15 -4.9941912711E+01 8.0000E+00 1.946E-05 4.863
16 -4.9941912691E+01 8.0000E+00 1.259E-05 4.785
17 -4.9941912688E+01 8.0000E+00 9.625E-06 4.602
18 -4.9941912717E+01 8.0000E+00 5.875E-06 4.505
19 -4.9941912717E+01 8.0000E+00 1.589E-06 4.619
20 -4.9941912716E+01 8.0000E+00 1.172E-06 4.574
21 -4.9941912703E+01 8.0000E+00 7.470E-07 4.412
1 -4.9866570408E+01 8.0000E+00 2.275E-01 19.253
2 -4.9871692332E+01 8.0000E+00 2.483E-01 6.240
3 -4.9931035176E+01 8.0000E+00 8.737E-02 3.726
4 -4.9938546535E+01 8.0000E+00 3.884E-02 3.758
5 -4.9935564525E+01 8.0000E+00 4.922E-02 3.564
6 -4.9941255002E+01 8.0000E+00 1.354E-02 3.622
7 -4.9941805662E+01 8.0000E+00 7.156E-03 3.557
8 -4.9941883346E+01 8.0000E+00 2.936E-03 3.551
9 -4.9941900213E+01 8.0000E+00 2.487E-03 3.263
10 -4.9941910449E+01 8.0000E+00 1.107E-03 3.380
11 -4.9941912266E+01 8.0000E+00 4.468E-04 3.414
12 -4.9941912597E+01 8.0000E+00 1.534E-04 3.523
13 -4.9941912657E+01 8.0000E+00 8.454E-05 3.111
14 -4.9941912680E+01 8.0000E+00 4.323E-05 3.022
15 -4.9941912712E+01 8.0000E+00 1.744E-05 3.127
16 -4.9941912692E+01 8.0000E+00 1.303E-05 3.053
17 -4.9941912689E+01 8.0000E+00 9.611E-06 3.000
18 -4.9941912713E+01 8.0000E+00 6.164E-06 2.883
19 -4.9941912717E+01 8.0000E+00 1.634E-06 2.917
20 -4.9941912715E+01 8.0000E+00 1.212E-06 2.979
21 -4.9941912704E+01 8.0000E+00 8.168E-07 2.789
Total number of SCF: 21
====================================================================
Energy and force calculation
====================================================================
Free energy per atom : -4.9941912703E+01 (Ha/atom)
Total free energy : -2.9965147622E+02 (Ha)
Band structure energy : -5.8070115209E+01 (Ha)
Exchange correlation energy : -4.9143438541E+01 (Ha)
Free energy per atom : -4.9941912704E+01 (Ha/atom)
Total free energy : -2.9965147623E+02 (Ha)
Band structure energy : -5.8070113340E+01 (Ha)
Exchange correlation energy : -4.9143437940E+01 (Ha)
Self and correction energy : -3.1654008538E+02 (Ha)
-Entropy*kb*T : -2.4357604837E-07 (Ha)
Fermi level : -1.9711942626E-01 (Ha)
RMS force : 5.4593647217E-03 (Ha/Bohr)
Maximum force : 7.9794015278E-03 (Ha/Bohr)
Time for force calculation : 0.348 (sec)
Maximum stress : 2.7162241740E+02 (Ha/Bohr**2)
Maximum stress equiv. to periodic : 2.6800506266E+05 (GPa)
Time for stress calculation : 0.348 (sec)
-Entropy*kb*T : -2.4357613067E-07 (Ha)
Fermi level : -1.9711941769E-01 (Ha)
Net magnetization : 7.9999800206E+00 (Bohr magneton)
RMS force : 5.4593856178E-03 (Ha/Bohr)
Maximum force : 7.9799382965E-03 (Ha/Bohr)
Time for force calculation : 0.312 (sec)
Maximum stress : 8.7178027121E-01 (Ha/Bohr)
Maximum stress equiv. to periodic : 3.3821423061E+00 (GPa)
Time for stress calculation : 0.427 (sec)
***************************************************************************
Timing info
***************************************************************************
Total walltime : 129.502 sec
Total walltime : 92.819 sec
___________________________________________________________________________

***************************************************************************
Expand Down
Loading
Loading