diff --git a/ChangeLog b/ChangeLog index 11d553d3..955af5c5 100644 --- a/ChangeLog +++ b/ChangeLog @@ -2,6 +2,12 @@ -Date -Name -changes +-------------- +September 04, 2025 +Name: Sayan Bhowmik +Changes: (src/md.c) +1. Bug fix in minimum pairwise distance calculation - include PBC/MIC. + -------------- September 03, 2025 Name: Tian Tian diff --git a/src/initialization.c b/src/initialization.c index 88ed6b6c..20ccb381 100644 --- a/src/initialization.c +++ b/src/initialization.c @@ -3722,7 +3722,7 @@ void write_output_init(SPARC_OBJ *pSPARC) { } fprintf(output_fp,"***************************************************************************\n"); - fprintf(output_fp,"* SPARC (version September 03, 2025) *\n"); + fprintf(output_fp,"* SPARC (version September 04, 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); diff --git a/src/md.c b/src/md.c index 0707fdc0..fc4e1445 100644 --- a/src/md.c +++ b/src/md.c @@ -13,6 +13,13 @@ #include #include +#ifdef USE_MKL + #include +#else + #include + #include +#endif + #include "md.h" #include "isddft.h" #include "orbitalElecDensInit.h" @@ -2135,11 +2142,68 @@ void MD_QOI(SPARC_OBJ *pSPARC, double *avgvel, double *maxvel, double *mindis) { //*avgdis = pow((pSPARC->range_x * pSPARC->range_y * pSPARC->range_z)/pSPARC->n_atom,1/3.0); int atm2, ityp2, cc2, pairIndex; cc = 0; + + // Store the cell as a 3x3 lattice vector + double *lattice = (double *)calloc(sizeof(double), 9); + double cell[3] = {pSPARC->range_x, pSPARC->range_y, pSPARC->range_z}; + double dr[3]; + int row, col; + for (row = 0; row < 3; row++) { + lattice[row*3] = pSPARC->LatUVec[row*3] * cell[row]; + lattice[row*3 + 1] = pSPARC->LatUVec[row*3 + 1] * cell[row]; + lattice[row*3 + 2] = pSPARC->LatUVec[row*3 + 2] * cell[row]; + } + + // Calculate the inverse of lattice-vector + double LV_inv[9], U[9], S[3], VT[9], superb[2], temp_mat[9], LatVec[9]; + double S_inv[9] = {0}; + int m = 3; + + for (int JJ = 0; JJ < 9; JJ++) { + LatVec[JJ] = lattice[JJ]; + } + + /* Calculating pinv(LatVec): This is necessary because for extremely skewed LV, the LV maybe close to singular */ + // Compute SVD of LatVec(LV) first: LV = U * S * V^T + LAPACKE_dgesvd(LAPACK_ROW_MAJOR, 'A', 'A', m, m, LatVec, m, S, U, m, VT, m, superb); + + // First compute S^{-1} + for (int i = 0; i < 3; i++) { + if (S[i] > 1e-12) S_inv[i * 3 + i] = 1.0/S[i]; + } + + // Compute LV^{-1} = V * S^{-1} * U^T + cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, m, m, m, 1.0, VT, m, S_inv, m, 0.0, temp_mat, m); + cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, m, m, m, 1.0, temp_mat, m, U, m, 0.0, LV_inv, m); + for(ityp = 0; ityp < pSPARC->Ntypes; ityp++){ mindis[ityp] = 1000000000.0; for(atm = 0; atm < pSPARC->nAtomv[ityp] - 1; atm++){ for(atm2 = atm + 1; atm2 < pSPARC->nAtomv[ityp]; atm2++){ - temp = fabs(sqrt(pow(pSPARC->atom_pos[3 * (atm + cc)] - pSPARC->atom_pos[3 * (atm2 + cc)],2.0) + pow(pSPARC->atom_pos[3 * (atm + cc) + 1] - pSPARC->atom_pos[3 * (atm2 + cc) + 1],2.0) + pow(pSPARC->atom_pos[3 * (atm + cc) + 2] - pSPARC->atom_pos[3 * (atm2 + cc) + 2],2.0) )); + // temp = fabs(sqrt(pow(pSPARC->atom_pos[3 * (atm + cc)] - pSPARC->atom_pos[3 * (atm2 + cc)],2.0) + pow(pSPARC->atom_pos[3 * (atm + cc) + 1] - pSPARC->atom_pos[3 * (atm2 + cc) + 1],2.0) + pow(pSPARC->atom_pos[3 * (atm + cc) + 2] - pSPARC->atom_pos[3 * (atm2 + cc) + 2],2.0) )); + + // Vector connecting atoms atm and atm2 + dr[0] = pSPARC->atom_pos[3 * (atm + cc)] - pSPARC->atom_pos[3 * (atm2 + cc)]; + dr[1] = pSPARC->atom_pos[3 * (atm + cc) + 1] - pSPARC->atom_pos[3 * (atm2 + cc) + 1]; + dr[2] = pSPARC->atom_pos[3 * (atm + cc) + 2] - pSPARC->atom_pos[3 * (atm2 + cc) + 2]; + + double frac[3], dr_pbc[3]; + + // Minimum Image Convention (MIC) in fractional coordinates + // frac = LV^{-1} * dr + cblas_dgemv(CblasRowMajor, CblasNoTrans, m, m, 1.0, LV_inv, m, dr, 1, 0.0, frac, 1); + + // Wrap into [-0.5, 0.5) + frac[0] -= round(frac[0]); + frac[1] -= round(frac[1]); + frac[2] -= round(frac[2]); + + // dr_pbc = lattice * frac + cblas_dgemv(CblasRowMajor, CblasNoTrans, m, m, 1.0, lattice, 3, frac, 1, 0.0, dr_pbc, 1); + + // Calculate distance + temp = sqrt(dr_pbc[0]*dr_pbc[0] + dr_pbc[1]*dr_pbc[1] + dr_pbc[2]*dr_pbc[2]); + if(temp < mindis[ityp]) mindis[ityp] = temp; } @@ -2155,7 +2219,30 @@ void MD_QOI(SPARC_OBJ *pSPARC, double *avgvel, double *maxvel, double *mindis) { 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) )); + // 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) )); + + // Vector connecting atoms atm and atm2 + dr[0] = pSPARC->atom_pos[3 * (atm + cc)] - pSPARC->atom_pos[3 * (atm2 + cc2)]; + dr[1] = pSPARC->atom_pos[3 * (atm + cc) + 1] - pSPARC->atom_pos[3 * (atm2 + cc2) + 1]; + dr[2] = pSPARC->atom_pos[3 * (atm + cc) + 2] - pSPARC->atom_pos[3 * (atm2 + cc2) + 2]; + + double frac[3], dr_pbc[3]; + + // Minimum Image Convention (MIC) in fractional coordinates + // frac = LV^{-1} * dr + cblas_dgemv(CblasRowMajor, CblasNoTrans, m, m, 1.0, LV_inv, m, dr, 1, 0.0, frac, 1); + + // Wrap into [-0.5, 0.5) + frac[0] -= round(frac[0]); + frac[1] -= round(frac[1]); + frac[2] -= round(frac[2]); + + // dr_pbc = lattice * frac + cblas_dgemv(CblasRowMajor, CblasNoTrans, m, m, 1.0, lattice, 3, frac, 1, 0.0, dr_pbc, 1); + + // Calculate distance + temp = sqrt(dr_pbc[0]*dr_pbc[0] + dr_pbc[1]*dr_pbc[1] + dr_pbc[2]*dr_pbc[2]); + if(temp < mindis[pairIndex]) mindis[pairIndex] = temp; } @@ -2165,6 +2252,7 @@ void MD_QOI(SPARC_OBJ *pSPARC, double *avgvel, double *maxvel, double *mindis) { } cc += pSPARC->nAtomv[ityp]; } + free(lattice); } /*