Skip to content

Commit 7436db8

Browse files
abhirajsharmaphanish-suryanarayana
abhirajsharma
authored and
phanish-suryanarayana
committed
Cyclix Stress bug fixed
1 parent 00dd1af commit 7436db8

42 files changed

Lines changed: 3682 additions & 3138 deletions

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎ChangeLog‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,14 @@
33
-Name
44
-changes
55
--------------
6+
September 25, 2025
7+
Name: Rajat kumar, Abhiraj Sharma
8+
Changes: src/cyclix/cyclix_stress.c, src/cyclix/cyclix_stress.h, src/stress.c, tests/cyclix
9+
1. Bug fixes in stress in cyclix due to more than 1 kpts/processor as well as spin calculation
10+
2. Bug fix in stress units in .out file as well as Ha/Bohr to GPa conversion
11+
3. Fixed reference results both standard as well as high accuracy in cyclix
12+
13+
614
September 04, 2025
715
Name: Sayan Bhowmik
816
Changes: (src/md.c)

‎src/cyclix/cyclix_stress.c‎

Lines changed: 20 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -853,17 +853,17 @@ void Compute_stress_tensor_kinetic_cyclix(SPARC_OBJ *pSPARC, double *dpsi, doubl
853853

854854
double temp_k = 0.0;
855855
for(n = 0; n < ncol; n++){
856-
double dpsi3_dpsi3 = 0.0;
857856
for (spinor = 0; spinor < Nspinor; spinor++) {
857+
double dpsi3_dpsi3 = 0.0;
858858
double *dpsi_ptr = dpsi + n * DMndsp + spinor * DMnd; // dpsi_1
859859
for(i = 0; i < DMnd; i++){
860860
dpsi3_dpsi3 += *(dpsi_ptr + i) * *(dpsi_ptr + i) * pSPARC->Intgwt_psi[i];
861861
}
862-
}
863-
double *occ = pSPARC->occ;
864-
if (pSPARC->spin_typ == 1) occ += spinor * Ns;
865-
double g_nk = occ[n + pSPARC->band_start_indx];
866-
temp_k += dpsi3_dpsi3 * g_nk;
862+
double *occ = pSPARC->occ;
863+
if (pSPARC->spin_typ == 1) occ += spinor * Ns;
864+
double g_nk = occ[n + pSPARC->band_start_indx];
865+
temp_k += dpsi3_dpsi3 * g_nk;
866+
}
867867
}
868868
*stress_k = - pSPARC->occfac * temp_k;
869869
}
@@ -882,7 +882,7 @@ void Compute_stress_tensor_nloc_by_integrals_cyclix(SPARC_OBJ *pSPARC, double *s
882882
IP_displ = pSPARC->IP_displ;
883883

884884
double *beta3_x3 = alpha + IP_displ[pSPARC->n_atom]*ncol*Nspinor;
885-
double stress_nl_ = 0;
885+
double stress_nl_ = 0.0;
886886
count = 0;
887887
for (spinor = 0; spinor < Nspinor; spinor++) {
888888
for (ityp = 0; ityp < pSPARC->Ntypes; ityp++) {
@@ -984,9 +984,9 @@ void Calculate_nonlocal_kinetic_stress_kpt_cyclix(SPARC_OBJ *pSPARC)
984984
beta3 = alpha_so2 + pSPARC->IP_displ_SOC[pSPARC->n_atom] * ncol * Nspinor * (Nk + kpt);
985985
Compute_Integral_Chi_XmRjp_beta_Dpsi_kpt_cyclix(pSPARC, pSPARC->Yorb_kpt, beta3, kpt, "SO2");
986986
}
987-
}
988987

989-
Compute_stress_tensor_kinetic_kpt_cyclix(pSPARC, pSPARC->Yorb_kpt, &stress_k);
988+
Compute_stress_tensor_kinetic_kpt_cyclix(pSPARC, pSPARC->Yorb_kpt, &stress_k, kpt);
989+
}
990990

991991
if (pSPARC->npNd > 1) {
992992
MPI_Allreduce(MPI_IN_PLACE, alpha, pSPARC->IP_displ[pSPARC->n_atom] * ncol * Nk * 2 * Nspinor, MPI_DOUBLE_COMPLEX, MPI_SUM, pSPARC->dmcomm);
@@ -1137,7 +1137,7 @@ void Compute_Integral_Chi_XmRjp_beta_Dpsi_kpt_cyclix(SPARC_OBJ *pSPARC, double _
11371137
/**
11381138
* @brief Calculate kinetic stress tensor
11391139
*/
1140-
void Compute_stress_tensor_kinetic_kpt_cyclix(SPARC_OBJ *pSPARC, double _Complex *dpsi, double *stress_k)
1140+
void Compute_stress_tensor_kinetic_kpt_cyclix(SPARC_OBJ *pSPARC, double _Complex *dpsi, double *stress_k, int kpt)
11411141
{
11421142
int ncol, DMnd, Nspinor, DMndsp, Nk, Ns;
11431143
ncol = pSPARC->Nband_bandcomm; // number of bands assigned
@@ -1147,26 +1147,22 @@ void Compute_stress_tensor_kinetic_kpt_cyclix(SPARC_OBJ *pSPARC, double _Complex
11471147
Nk = pSPARC->Nkpts_kptcomm;
11481148
Ns = pSPARC->Nstates;
11491149

1150-
int kpt, n, spinor, i;
1151-
double stress_k_ = 0;
1152-
for(kpt = 0; kpt < Nk; kpt++) {
1153-
double temp_k = 0.0;
1154-
for(n = 0; n < ncol; n++){
1150+
int n, spinor, i;
1151+
double temp_k = 0.0;
1152+
for(n = 0; n < ncol; n++){
1153+
for (spinor = 0; spinor < Nspinor; spinor++) {
11551154
double dpsi3_dpsi3 = 0.0;
1156-
for (spinor = 0; spinor < Nspinor; spinor++) {
1157-
double _Complex *dpsi_ptr = dpsi + n * DMndsp + spinor * DMnd; // dpsi_1
1158-
for(i = 0; i < DMnd; i++){
1159-
dpsi3_dpsi3 += (creal(*(dpsi_ptr + i)) * creal(*(dpsi_ptr + i)) + cimag(*(dpsi_ptr + i)) * cimag(*(dpsi_ptr + i))) * pSPARC->Intgwt_psi[i];
1160-
}
1155+
double _Complex *dpsi_ptr = dpsi + n * DMndsp + spinor * DMnd; // dpsi_1
1156+
for(i = 0; i < DMnd; i++){
1157+
dpsi3_dpsi3 += (creal(*(dpsi_ptr + i)) * creal(*(dpsi_ptr + i)) + cimag(*(dpsi_ptr + i)) * cimag(*(dpsi_ptr + i))) * pSPARC->Intgwt_psi[i];
11611158
}
11621159
double *occ = pSPARC->occ + kpt*Ns;
11631160
if (pSPARC->spin_typ == 1) occ += spinor * Nk * Ns;
11641161
double g_nk = occ[n + pSPARC->band_start_indx];
11651162
temp_k += dpsi3_dpsi3 * g_nk;
11661163
}
1167-
stress_k_ -= pSPARC->occfac * pSPARC->kptWts_loc[kpt] / pSPARC->Nkpts * temp_k;
11681164
}
1169-
*stress_k = stress_k_;
1165+
*stress_k -= pSPARC->occfac * pSPARC->kptWts_loc[kpt] / pSPARC->Nkpts * temp_k;
11701166
}
11711167

11721168

@@ -1175,7 +1171,7 @@ void Compute_stress_tensor_nloc_by_integrals_kpt_cyclix(SPARC_OBJ *pSPARC, doubl
11751171
int k, n, np, ldispl, ityp, iat, ncol, Ns;
11761172
int count, l, m, lmax, Nk, spinor, Nspinor;
11771173
int l_start, mexclude, ppl, *IP_displ;
1178-
double g_nk, kptwt, alpha_r, alpha_i, SJ, temp2_s, scaled_gamma_Jl = 0;
1174+
double g_nk, kptwt, alpha_r, alpha_i, SJ, temp2_s, scaled_gamma_Jl = 0.0;
11791175
ncol = pSPARC->Nband_bandcomm; // number of bands assigned
11801176
Ns = pSPARC->Nstates;
11811177
Nk = pSPARC->Nkpts_kptcomm;
@@ -1233,4 +1229,4 @@ void Compute_stress_tensor_nloc_by_integrals_kpt_cyclix(SPARC_OBJ *pSPARC, doubl
12331229
}
12341230
}
12351231
}
1236-
}
1232+
}

‎src/cyclix/include/cyclix_stress.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,7 @@ void Dpseudopot_cyclix_z(SPARC_OBJ *pSPARC, double *VJ, int FDn, int ishift_p, i
4545
* @brief Calculate kinetic stress tensor
4646
*/
4747
void Compute_stress_tensor_kinetic_cyclix(SPARC_OBJ *pSPARC, double *dpsi, double *stress_k);
48-
void Compute_stress_tensor_kinetic_kpt_cyclix(SPARC_OBJ *pSPARC, double _Complex *dpsi, double *stress_k);
48+
void Compute_stress_tensor_kinetic_kpt_cyclix(SPARC_OBJ *pSPARC, double _Complex *dpsi, double *stress_k, int kpt);
4949

5050
/**
5151
* @brief Calculate nonlocal stress components

‎src/electronicGroundState.c‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -284,7 +284,7 @@ void Calculate_electronicGroundState(SPARC_OBJ *pSPARC) {
284284
double cellsizes[3] = {pSPARC->range_x, pSPARC->range_y, pSPARC->range_z};
285285
int BCs[3] = {pSPARC->BCx, pSPARC->BCy, pSPARC->BCz};
286286
char stressUnit[16];
287-
double maxSGPa = convertStressToGPa(maxS, cellsizes, BCs, stressUnit);
287+
double maxSGPa = convertStressToGPa(pSPARC, maxS, cellsizes, BCs, stressUnit);
288288
fprintf(output_fp,"Maximum stress :%18.10E (%s)\n",maxS,stressUnit);
289289
fprintf(output_fp,"Maximum stress equiv. to periodic :%18.10E (GPa)\n",maxSGPa);
290290
}

‎src/include/stress.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -92,7 +92,7 @@ void PrintStress (SPARC_OBJ *pSPARC, double *stress, FILE *fp);
9292
* @param origUnit (OUTPUT) Unit for the original stress.
9393
* @return double Stress component in GPa.
9494
*/
95-
double convertStressToGPa(double Stress, double cellsizes[3], int BCs[3], char origUnit[16]);
95+
double convertStressToGPa(SPARC_OBJ *pSPARC, double Stress, double cellsizes[3], int BCs[3], char origUnit[16]);
9696

9797

9898
#endif // STRESS_H

‎src/initialization.c‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3722,7 +3722,7 @@ void write_output_init(SPARC_OBJ *pSPARC) {
37223722
}
37233723

37243724
fprintf(output_fp,"***************************************************************************\n");
3725-
fprintf(output_fp,"* SPARC (version September 04, 2025) *\n");
3725+
fprintf(output_fp,"* SPARC (version September 25, 2025) *\n");
37263726
fprintf(output_fp,"* Copyright (c) 2020 Material Physics & Mechanics Group, Georgia Tech *\n");
37273727
fprintf(output_fp,"* Distributed under GNU General Public License 3 (GPL) *\n");
37283728
fprintf(output_fp,"* Start time: %s *\n",c_time_str);

‎src/makefile‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -159,4 +159,4 @@ clean:
159159
rm -f $(OBJSC) $(LIBBASE)
160160
rm -f socket/*.o
161161
test: ../tests/SPARC_testing_script.py
162-
cd ../tests; python SPARC_testing_script.py
162+
cd ../tests; python SPARC_testing_script.py

‎src/stress.c‎

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -2113,22 +2113,25 @@ void PrintStress (SPARC_OBJ *pSPARC, double *stress, FILE *fp) {
21132113
* @param origUnit (OUTPUT) Unit for the original stress.
21142114
* @return double Stress component in GPa.
21152115
*/
2116-
double convertStressToGPa(double Stress, double cellsizes[3], int BCs[3], char origUnit[16]) {
2116+
double convertStressToGPa(SPARC_OBJ *pSPARC, double Stress, double cellsizes[3], int BCs[3], char origUnit[16]) {
21172117
double scale = 1.0; // scale the stress by the dimensions of the dirichlet directions
21182118
int nperiods = 3;
21192119
if (BCs[0] == 1) { scale *= cellsizes[0]; nperiods--; }
21202120
if (BCs[1] == 1) { scale *= cellsizes[1]; nperiods--; }
21212121
if (BCs[2] == 1) { scale *= cellsizes[2]; nperiods--; }
21222122

21232123
// scale and then convert to GPa
2124-
double StressGPa = Stress / scale * CONST_HA_BOHR3_GPA;
2124+
double StressGPa;
2125+
if (pSPARC->CyclixFlag)
2126+
StressGPa = Stress*CONST_HA_BOHR3_GPA/(M_PI * ( pow(pSPARC->xout,2.0) - pow(pSPARC->xin,2.0) ) );
2127+
else
2128+
StressGPa = Stress / scale * CONST_HA_BOHR3_GPA;
21252129

21262130
// find the original unit
2127-
if (nperiods == 1)
2131+
if (nperiods == 1 || pSPARC->CyclixFlag)
21282132
snprintf(origUnit, 16, "Ha/Bohr");
21292133
else
21302134
snprintf(origUnit, 16, "Ha/Bohr**%d", nperiods);
2131-
21322135
return StressGPa;
21332136
}
21342137

‎tests/SPARC_testing_script.py‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -12,12 +12,12 @@
1212
import math
1313

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

2222

2323

‎tests/cyclix/FeCl2_cyclix_spin/high_accuracy/FeCl2_cyclix_spin.refout‎

Lines changed: 45 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,8 @@
11
***************************************************************************
2-
* SPARC (version June 24, 2024) *
2+
* SPARC (version September 04, 2025) *
33
* Copyright (c) 2020 Material Physics & Mechanics Group, Georgia Tech *
44
* Distributed under GNU General Public License 3 (GPL) *
5-
* Start time: Mon Jun 24 15:30:15 2024 *
5+
* Start time: Wed Sep 24 20:42:08 2025 *
66
***************************************************************************
77
Input parameters
88
***************************************************************************
@@ -17,6 +17,7 @@ SPIN_TYP: 1
1717
ELEC_TEMP_TYPE: Fermi-Dirac
1818
SMEARING: 0.001
1919
EXCHANGE_CORRELATION: GGA_PBE
20+
HUBBARD_FLAG: 0
2021
NSTATES: 41
2122
CHEB_DEGREE: 100
2223
CHEFSI_BOUND_FLAG: 0
@@ -47,7 +48,7 @@ PRINT_ATOMS: 1
4748
PRINT_EIGEN: 0
4849
PRINT_DENSITY: 0
4950
PRINT_ENERGY_DENSITY: 0
50-
OUTPUT_FILE: FeCl2_cyclix_spin/temp_run/FeCl2_cyclix_spin
51+
OUTPUT_FILE: cyclix/FeCl2_cyclix_spin/temp_run/FeCl2_cyclix_spin
5152
***************************************************************************
5253
Cell
5354
***************************************************************************
@@ -58,18 +59,18 @@ Density: 1.2141297401E-01 (amu/Bohr^3), 1.3605383749E+00 (g/cc)
5859
***************************************************************************
5960
NP_SPIN_PARAL: 2
6061
NP_KPOINT_PARAL: 1
61-
NP_BAND_PARAL: 6
62-
NP_DOMAIN_PARAL: 4 1 2
63-
NP_DOMAIN_PHI_PARAL: 12 2 4
62+
NP_BAND_PARAL: 7
63+
NP_DOMAIN_PARAL: 17 1 1
64+
NP_DOMAIN_PHI_PARAL: 16 3 5
6465
EIG_SERIAL_MAXNS: 1500
6566
***************************************************************************
6667
Initialization
6768
***************************************************************************
68-
Number of processors : 96
69+
Number of processors : 240
6970
Mesh spacing in x-direction : 0.120234 (Bohr)
7071
Mesh spacing in y-direction : 0.00296377 (Bohr)
7172
Mesh spacing in z-direction : 0.119703 (Bohr)
72-
Output printed to : FeCl2_cyclix_spin/temp_run/FeCl2_cyclix_spin.out
73+
Output printed to : cyclix/FeCl2_cyclix_spin/temp_run/FeCl2_cyclix_spin.out
7374
Total number of atom types : 2
7475
Total number of atoms : 6
7576
Total number of electrons : 60
@@ -84,53 +85,54 @@ Atomic mass : 35.4515
8485
Pseudocharge radii of atom type 2 : 6.61 0.15 6.58 (x, y, z dir)
8586
Number of atoms of type 2 : 4
8687
Estimated total memory usage : 5.82 GB
87-
Estimated memory per processor : 62.09 MB
88+
Estimated memory per processor : 24.84 MB
8889
========================================================================================
8990
Self Consistent Field (SCF#1)
9091
========================================================================================
9192
Iteration Free Energy (Ha/atom) Magnetization SCF Error Timing (sec)
92-
1 -4.9866567407E+01 8.0000E+00 2.275E-01 18.758
93-
2 -4.9871693480E+01 8.0000E+00 2.483E-01 6.083
94-
3 -4.9931035419E+01 8.0000E+00 8.737E-02 5.879
95-
4 -4.9938547782E+01 8.0000E+00 3.884E-02 5.766
96-
5 -4.9935561172E+01 8.0000E+00 4.924E-02 5.708
97-
6 -4.9941254720E+01 8.0000E+00 1.354E-02 5.598
98-
7 -4.9941805651E+01 8.0000E+00 7.158E-03 5.669
99-
8 -4.9941883346E+01 8.0000E+00 2.936E-03 5.329
100-
9 -4.9941900215E+01 8.0000E+00 2.487E-03 5.001
101-
10 -4.9941910449E+01 8.0000E+00 1.106E-03 5.320
102-
11 -4.9941912265E+01 8.0000E+00 4.472E-04 8.508
103-
12 -4.9941912597E+01 8.0000E+00 1.529E-04 5.152
104-
13 -4.9941912656E+01 8.0000E+00 8.555E-05 4.903
105-
14 -4.9941912686E+01 8.0000E+00 4.390E-05 4.870
106-
15 -4.9941912711E+01 8.0000E+00 1.946E-05 4.863
107-
16 -4.9941912691E+01 8.0000E+00 1.259E-05 4.785
108-
17 -4.9941912688E+01 8.0000E+00 9.625E-06 4.602
109-
18 -4.9941912717E+01 8.0000E+00 5.875E-06 4.505
110-
19 -4.9941912717E+01 8.0000E+00 1.589E-06 4.619
111-
20 -4.9941912716E+01 8.0000E+00 1.172E-06 4.574
112-
21 -4.9941912703E+01 8.0000E+00 7.470E-07 4.412
93+
1 -4.9866570408E+01 8.0000E+00 2.275E-01 19.253
94+
2 -4.9871692332E+01 8.0000E+00 2.483E-01 6.240
95+
3 -4.9931035176E+01 8.0000E+00 8.737E-02 3.726
96+
4 -4.9938546535E+01 8.0000E+00 3.884E-02 3.758
97+
5 -4.9935564525E+01 8.0000E+00 4.922E-02 3.564
98+
6 -4.9941255002E+01 8.0000E+00 1.354E-02 3.622
99+
7 -4.9941805662E+01 8.0000E+00 7.156E-03 3.557
100+
8 -4.9941883346E+01 8.0000E+00 2.936E-03 3.551
101+
9 -4.9941900213E+01 8.0000E+00 2.487E-03 3.263
102+
10 -4.9941910449E+01 8.0000E+00 1.107E-03 3.380
103+
11 -4.9941912266E+01 8.0000E+00 4.468E-04 3.414
104+
12 -4.9941912597E+01 8.0000E+00 1.534E-04 3.523
105+
13 -4.9941912657E+01 8.0000E+00 8.454E-05 3.111
106+
14 -4.9941912680E+01 8.0000E+00 4.323E-05 3.022
107+
15 -4.9941912712E+01 8.0000E+00 1.744E-05 3.127
108+
16 -4.9941912692E+01 8.0000E+00 1.303E-05 3.053
109+
17 -4.9941912689E+01 8.0000E+00 9.611E-06 3.000
110+
18 -4.9941912713E+01 8.0000E+00 6.164E-06 2.883
111+
19 -4.9941912717E+01 8.0000E+00 1.634E-06 2.917
112+
20 -4.9941912715E+01 8.0000E+00 1.212E-06 2.979
113+
21 -4.9941912704E+01 8.0000E+00 8.168E-07 2.789
113114
Total number of SCF: 21
114115
====================================================================
115116
Energy and force calculation
116117
====================================================================
117-
Free energy per atom : -4.9941912703E+01 (Ha/atom)
118-
Total free energy : -2.9965147622E+02 (Ha)
119-
Band structure energy : -5.8070115209E+01 (Ha)
120-
Exchange correlation energy : -4.9143438541E+01 (Ha)
118+
Free energy per atom : -4.9941912704E+01 (Ha/atom)
119+
Total free energy : -2.9965147623E+02 (Ha)
120+
Band structure energy : -5.8070113340E+01 (Ha)
121+
Exchange correlation energy : -4.9143437940E+01 (Ha)
121122
Self and correction energy : -3.1654008538E+02 (Ha)
122-
-Entropy*kb*T : -2.4357604837E-07 (Ha)
123-
Fermi level : -1.9711942626E-01 (Ha)
124-
RMS force : 5.4593647217E-03 (Ha/Bohr)
125-
Maximum force : 7.9794015278E-03 (Ha/Bohr)
126-
Time for force calculation : 0.348 (sec)
127-
Maximum stress : 2.7162241740E+02 (Ha/Bohr**2)
128-
Maximum stress equiv. to periodic : 2.6800506266E+05 (GPa)
129-
Time for stress calculation : 0.348 (sec)
123+
-Entropy*kb*T : -2.4357613067E-07 (Ha)
124+
Fermi level : -1.9711941769E-01 (Ha)
125+
Net magnetization : 7.9999800206E+00 (Bohr magneton)
126+
RMS force : 5.4593856178E-03 (Ha/Bohr)
127+
Maximum force : 7.9799382965E-03 (Ha/Bohr)
128+
Time for force calculation : 0.312 (sec)
129+
Maximum stress : 8.7178027121E-01 (Ha/Bohr)
130+
Maximum stress equiv. to periodic : 3.3821423061E+00 (GPa)
131+
Time for stress calculation : 0.427 (sec)
130132
***************************************************************************
131133
Timing info
132134
***************************************************************************
133-
Total walltime : 129.502 sec
135+
Total walltime : 92.819 sec
134136
___________________________________________________________________________
135137

136138
***************************************************************************

0 commit comments

Comments
 (0)