Skip to content

Commit

Permalink
Fixing bug in potential
Browse files Browse the repository at this point in the history
  • Loading branch information
ssahoo41 committed Jan 22, 2024
1 parent ae0e5e5 commit 0242c6d
Show file tree
Hide file tree
Showing 6 changed files with 55 additions and 170 deletions.
94 changes: 47 additions & 47 deletions src/convolutional_gga/convolutionalGGAMultipoleXC.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,9 +213,9 @@ void Calculate_Vx_GGA_CONV_PBE_MULTIPOLE(SPARC_OBJ *pSPARC, XCCST_OBJ *xc_cst, d
dfxdg = dfxdss * dssdg;
Dxdgrho[i] = 0.5 * ex_lsd * rho_updn * dfxdg;// second part of the derivative
ex = ex_gga * rho_updn;
//Df_alpha = -xc_cst->kappa * pow(divss, alpha[i])*(-log(1/divss) + (divss/alpha[i]));
//Dalpha_qp = -((4.0-0.75) * m)/(2.0 + exp(- m * (feat_qp_monopole[i] - n)) + exp(m * (feat_qp_monopole[i] - n)));
//Dfeat_qp_mp[i] = ex_lsd * 2 * rho_updn * Df_alpha * Dalpha_qp;
Df_alpha = -xc_cst->kappa * pow(divss, alpha[i])*(-log(1/divss) + ((xc_cst->mu_divkappa * ss * divss)/alpha[i])); // fixing bug (deleted the fix for sometime)
Dalpha_qp = -((4.0-0.75) * m)/(2.0 + exp(- m * (feat_qp_monopole[i] - n)) + exp(m * (feat_qp_monopole[i] - n)));
Dfeat_qp_mp[i] = ex_lsd * 2 * rho_updn * Df_alpha * Dalpha_qp;

// If non spin-polarized, treat spin down contribution now, similar to spin up
ex = ex * 2.0;
Expand Down Expand Up @@ -263,38 +263,38 @@ void Calculate_Vx_GGA_CONV_PBE_MULTIPOLE(SPARC_OBJ *pSPARC, XCCST_OBJ *xc_cst, d
free(pbeX_potential);
#endif

// // Extra potential contribution from HSMP part saved to a different file
// #ifdef DEBUG
// if (rank == 0) printf("Saving extra potential contribution from HSMP.\n");
// #endif
// double *extra_potential, *global_extra_potential;
// extra_potential = (double *)malloc(DMnd * sizeof(double));
// for (i = 0; i < DMnd; i++){
// extra_potential[i] = 0;
// }
// #ifdef DEBUG
// global_extra_potential = (double *)malloc(pSPARC->Nd * sizeof(double));
// #endif

// gather_distributed_vector(Dfeat_qp_mp, pSPARC->DMVertices, global_Df_featmp, gridsizes, pSPARC->dmcomm_phi, 1);
// MPI_Bcast(global_Df_featmp, pSPARC->Nd, MPI_DOUBLE, 0, MPI_COMM_WORLD);
// Conv_feat_vectors_dir(pSPARC, &mp, pSPARC->DMVertices, 1, global_Df_featmp, DDfeat_qp_mp, "000", pSPARC->dmcomm_phi);
// // monopole contribution
// for (i = 0; i < DMnd; i++){
// XPotential[i] += DDfeat_qp_mp[i];
// extra_potential[i] += DDfeat_qp_mp[i];
// }

// #ifdef DEBUG
// gather_distributed_vector(extra_potential, pSPARC->DMVertices, global_extra_potential, gridsizes, pSPARC->dmcomm_phi, 1);
// char extraPotentialFilename[128];
// if (rank == 0){
// snprintf(extraPotentialFilename, 128, "extra_X_potential.csv");
// writeMatToFile(extraPotentialFilename, global_extra_potential, pSPARC->Nx, pSPARC->Ny, pSPARC->Nz);
// }
// #endif
// Extra potential contribution from HSMP part saved to a different file
#ifdef DEBUG
if (rank == 0) printf("Saving extra potential contribution from HSMP.\n");
#endif
double *extra_potential, *global_extra_potential;
extra_potential = (double *)malloc(DMnd * sizeof(double));
for (i = 0; i < DMnd; i++){
extra_potential[i] = 0;
}
#ifdef DEBUG
global_extra_potential = (double *)malloc(pSPARC->Nd * sizeof(double));
#endif

gather_distributed_vector(Dfeat_qp_mp, pSPARC->DMVertices, global_Df_featmp, gridsizes, pSPARC->dmcomm_phi, 1);
MPI_Bcast(global_Df_featmp, pSPARC->Nd, MPI_DOUBLE, 0, MPI_COMM_WORLD);
Conv_feat_vectors_dir(pSPARC, &mp, pSPARC->DMVertices, 1, global_Df_featmp, DDfeat_qp_mp, "000", pSPARC->dmcomm_phi);
// monopole contribution
for (i = 0; i < DMnd; i++){
XPotential[i] += DDfeat_qp_mp[i];
extra_potential[i] += DDfeat_qp_mp[i];
}

#ifdef DEBUG
gather_distributed_vector(extra_potential, pSPARC->DMVertices, global_extra_potential, gridsizes, pSPARC->dmcomm_phi, 1);
char extraPotentialFilename[128];
if (rank == 0){
snprintf(extraPotentialFilename, 128, "extra_X_potential.csv");
writeMatToFile(extraPotentialFilename, global_extra_potential, pSPARC->Nx, pSPARC->Ny, pSPARC->Nz);
}
#endif
//Deallocate memory
// free(extra_potential); free(global_extra_potential);
free(extra_potential); free(global_extra_potential);
free(Drho_x); free(Drho_y); free(Drho_z);
free(global_rho);
free(DDrho_x); free(DDrho_y); free(DDrho_z);
Expand Down Expand Up @@ -489,9 +489,9 @@ void Calculate_Vx_GSGA_CONV_PBE_MULTIPOLE(SPARC_OBJ *pSPARC, XCCST_OBJ *xc_cst,
dfxdg = dfxdss * dssdg;
Dxdgrho[spn_i*DMnd + i] = ex_lsd * rho_updn * dfxdg; // spin up and spin down for second part of the derivative
ex += ex_gga * rho_updn;
// Df_alpha = -xc_cst->kappa * pow(divss, alpha[spn_i*DMnd + i])*(-log(1/divss) + (divss/alpha[spn_i*DMnd + i]));
// Dalpha_qp = - ((4.0 - 0.75) * m)/(2.0 + exp(-m * (feat_qp_monopole[spn_i*DMnd + i]-n)) + exp(m*(feat_qp_monopole[spn_i*DMnd + i]-n)));
// Dfeat_qp_mp[spn_i*DMnd+i] = ex_lsd * rho_updn * Df_alpha * Dalpha_qp;
Df_alpha = -xc_cst->kappa * pow(divss, alpha[spn_i*DMnd + i])*(-log(1/divss) + ((xc_cst->mu_divkappa * ss * divss)/alpha[spn_i*DMnd + i]));
Dalpha_qp = - ((4.0 - 0.75) * m)/(2.0 + exp(-m * (feat_qp_monopole[spn_i*DMnd + i]-n)) + exp(m*(feat_qp_monopole[spn_i*DMnd + i]-n)));
Dfeat_qp_mp[spn_i*DMnd+i] = ex_lsd * rho_updn * Df_alpha * Dalpha_qp;
}
e_x[i] = ex * rhotot_inv;
}
Expand Down Expand Up @@ -522,18 +522,18 @@ void Calculate_Vx_GSGA_CONV_PBE_MULTIPOLE(SPARC_OBJ *pSPARC, XCCST_OBJ *xc_cst,
XPotential[DMnd + i] += - DDrho_x[DMnd + i] - DDrho_y[DMnd + i] - DDrho_z[DMnd + i];
}

// // add the contribution from the monopole
// for (int i = 0; i < 2; i++){
// gather_distributed_vector(Dfeat_qp_mp + i*DMnd, pSPARC->DMVertices, global_Df_featmp + i*Nd, gridsizes, pSPARC->dmcomm_phi, 1);
// }
// add the contribution from the monopole
for (int i = 0; i < 2; i++){
gather_distributed_vector(Dfeat_qp_mp + i*DMnd, pSPARC->DMVertices, global_Df_featmp + i*Nd, gridsizes, pSPARC->dmcomm_phi, 1);
}

// MPI_Bcast(global_Df_featmp, 2*Nd, MPI_DOUBLE, 0, pSPARC->dmcomm_phi);
// Conv_feat_vectors_dir(pSPARC, &mp, pSPARC->DMVertices, 2, global_Df_featmp, DDfeat_qp_mp, "000", pSPARC->dmcomm_phi);
MPI_Bcast(global_Df_featmp, 2*Nd, MPI_DOUBLE, 0, pSPARC->dmcomm_phi);
Conv_feat_vectors_dir(pSPARC, &mp, pSPARC->DMVertices, 2, global_Df_featmp, DDfeat_qp_mp, "000", pSPARC->dmcomm_phi);

// for(i=0; i < DMnd; i++){
// XPotential[i] += DDfeat_qp_mp[i];
// XPotential[DMnd + i] += DDfeat_qp_mp[DMnd + i];
// }
for(i=0; i < DMnd; i++){
XPotential[i] += DDfeat_qp_mp[i];
XPotential[DMnd + i] += DDfeat_qp_mp[DMnd + i];
}


// Deallocate memory
Expand Down
7 changes: 0 additions & 7 deletions src/exchangeCorrelation.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -443,10 +443,6 @@ void Calculate_Vxc_GGA_PBE(SPARC_OBJ *pSPARC, XCCST_OBJ *xc_cst, double *rho) {
if (pSPARC->dmcomm_phi == MPI_COMM_NULL) {
return;
}

//double beta, kappa_pbe, mu, rsfac, kappa, mu_divkappa_pbe, mu_divkappa;
//double ec0_aa, ec0_a1, ec0_b1, ec0_b2, ec0_b3, ec0_b4;
//double pi, piinv, third, twom1_3, sixpi2_1_3, sixpi2m1_3, threefourth_divpi, gamma, gamma_inv, coeff_tt, sq_rsfac, sq_rsfac_inv;
double rho_updn, rho_updnm1_3, rhom1_3, rhotot_inv, rhotmo6, rhoto6, rhomot, ex_lsd, rho_inv, coeffss, ss;
double divss, dfxdss, fx, ex_gga, dssdn, dfxdn, dssdg, dfxdg, exc, rs, sqr_rs, rsm1_2;
double ec0_q0, ec0_q1, ec0_q1p, ec0_den, ec0_log, ecrs0, decrs0_drs, ecrs, decrs_drs;
Expand Down Expand Up @@ -486,13 +482,10 @@ void Calculate_Vxc_GGA_PBE(SPARC_OBJ *pSPARC, XCCST_OBJ *xc_cst, double *rho) {
}

// Compute exchange and correlation
//phi_zeta = 1.0;
//phip_zeta = 0.0;
phi_zeta_inv = 1.0;
//phi_logder = 0.0;
phi3_zeta = 1.0;
gamphi3inv = xc_cst->gamma_inv;
//phipp_zeta = (-2.0/9.0) * alpha_zeta * alpha_zeta;

for(i = 0; i < DMnd; i++){
rho_updn = rho[i]/2.0;
Expand Down
1 change: 0 additions & 1 deletion src/multipole_features/MCSH.c
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@

/**
* @file MCSH.c
* @brief This file contains the declaration of MCSH functions.
Expand Down
27 changes: 7 additions & 20 deletions src/multipole_features/MCSHHelper.c
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
/**
* @file MCSHHelper.h
* @file MCSHHelper.c
* @brief This file contains the declaration of MCSH functions.
*
* @author Sushree Jagriti Sahoo <[email protected]>
Expand Down Expand Up @@ -335,8 +335,7 @@ void getRArray(const double *x, const double *y, const double *z, double *result
void polyArray(const double *x, const int powX, const double a, double *result, const int size)
{
int i = 0;
for (i = 0; i < size; i++)
{
for (i = 0; i < size; i++){
result[i] = a * pow(x[i], powX);
}
}
Expand All @@ -347,8 +346,7 @@ void polyArray(const double *x, const int powX, const double a, double *result,
void polyArray_add(const double *x, const int powX, const double a, double *result, const int size)
{
int i = 0;
for (i = 0; i < size; i++)
{
for (i = 0; i < size; i++){
result[i] += a * pow(x[i], powX);
}
}
Expand All @@ -359,8 +357,7 @@ void polyArray_add(const double *x, const int powX, const double a, double *resu
void polyXYZArray(const double *x, const double *y, const double *z, const int powX, const int powY, const int powZ, const double a, double *result, const int size)
{
int i = 0;
for (i = 0; i < size; i++)
{
for (i = 0; i < size; i++){
result[i] = a * pow(x[i], powX) * pow(y[i], powY) * pow(z[i], powZ) ;
}
}
Expand All @@ -371,8 +368,7 @@ void polyXYZArray(const double *x, const double *y, const double *z, const int p
void polyXYZArray_add(const double *x, const double *y, const double *z, const int powX, const int powY, const int powZ, const double a, double *result, const int size)
{
int i = 0;
for (i = 0; i < size; i++)
{
for (i = 0; i < size; i++){
result[i] += a * pow(x[i], powX) * pow(y[i], powY) * pow(z[i], powZ) ;
}
}
Expand All @@ -386,15 +382,13 @@ void applyU(double *X, double *Y, double*Z, const double *U, const int size)
double *matMutResult = malloc( 3 * size * sizeof(double));

int i;
for ( i = 0; i < size; i++ )
{
for ( i = 0; i < size; i++ ){
combinedArr[ i * 3 ] = X[i];
combinedArr[ i * 3 + 1 ] = Y[i];
combinedArr[ i * 3 + 2 ] = Z[i];
}

for ( i = 0; i < size; i++ )
{
for ( i = 0; i < size; i++ ){
X[i] = matMutResult[ i * 3 ];
Y[i] = matMutResult[ i * 3 + 1 ];
Z[i] = matMutResult[ i * 3 + 2 ];
Expand Down Expand Up @@ -1063,7 +1057,6 @@ void getMainParameter_RadialLegendre(MULTIPOLE_OBJ *mp, int* LegendreOrderList,
{
LegendreOrderList[index] = j;
lList[index] = getCurrentLNumber(i);
//groupList[index] = getCurrentGroupNumber(i);
index++;
}
}
Expand All @@ -1074,15 +1067,9 @@ void getMainParameter_RadialLegendre(MULTIPOLE_OBJ *mp, int* LegendreOrderList,
*/
int getDescriptorListLength_RadialLegendre(MULTIPOLE_OBJ *mp)
{

//int numGroup = getNumGroup(maxMCSHOrder);
// printf("\nnumber groups:%d \n", numGroup);
int numMCSH = mp->MCSHMaxOrder + 1;
int numLegendre = mp->MCSHMaxRadialOrder + 1;
// printf("\nnumber r cutoff:%d \n", numRCutoff);

return numMCSH * numLegendre;

}

/**
Expand Down
89 changes: 0 additions & 89 deletions src/multipole_features/MCSHMainCalc.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -226,95 +226,6 @@ void Multipole_Initialize(SPARC_OBJ *pSPARC, MULTIPOLE_OBJ *mp) {
}
}

/**
* @brief function to collect distributed rho vector into a global vector
*/
void CollectElecDens(SPARC_OBJ *pSPARC, double *global_rho, double *global_psdchrgdens) {
if (pSPARC->dmcomm_phi == MPI_COMM_NULL) return;
int nproc_dmcomm_phi, rank_dmcomm_phi, DMnd, i, j, k, index;
MPI_Comm_size(pSPARC->dmcomm_phi, &nproc_dmcomm_phi);
MPI_Comm_rank(pSPARC->dmcomm_phi, &rank_dmcomm_phi);

int Nd = pSPARC->Nd;
DMnd = pSPARC->Nd_d;

double *rho, *b;
rho = NULL;
b = NULL;
if (nproc_dmcomm_phi > 1) { // if there's more than one process, need to collect rho first
// use DD2DD to collect distributed data
int gridsizes[3], sdims[3], rdims[3], rDMVert[6];
MPI_Comm recv_comm;
if (rank_dmcomm_phi) {
recv_comm = MPI_COMM_NULL;
} else {
int dims[3] = {1,1,1}, periods[3] = {1,1,1};
// create a cartesian topology on one process (rank 0)
MPI_Cart_create(MPI_COMM_SELF, 3, dims, periods, 0, &recv_comm);
}
D2D_OBJ d2d_sender, d2d_recvr;
gridsizes[0] = pSPARC->Nx;
gridsizes[1] = pSPARC->Ny;
gridsizes[2] = pSPARC->Nz;
sdims[0] = pSPARC->npNdx_phi;
sdims[1] = pSPARC->npNdy_phi;
sdims[2] = pSPARC->npNdz_phi;
rdims[0] = rdims[1] = rdims[2] = 1;
rDMVert[0] = 0; rDMVert[1] = pSPARC->Nx-1;
rDMVert[2] = 0; rDMVert[3] = pSPARC->Ny-1;
rDMVert[4] = 0; rDMVert[5] = pSPARC->Nz-1;

// set up D2D targets
Set_D2D_Target(&d2d_sender, &d2d_recvr, gridsizes, pSPARC->DMVertices, rDMVert, pSPARC->dmcomm_phi,
sdims, recv_comm, rdims, pSPARC->dmcomm_phi);
if (rank_dmcomm_phi == 0) {
int n_rho = pSPARC->Nspden/2*2+1;
rho = (double*)malloc(pSPARC->Nd * n_rho * sizeof(double)); // allocating memory in only one processor
b = (double*)malloc(pSPARC->Nd * sizeof(double));
}

// send rho

D2D(&d2d_sender, &d2d_recvr, gridsizes, pSPARC->DMVertices, pSPARC->electronDens, rDMVert,
rho, pSPARC->dmcomm_phi, sdims, recv_comm, rdims, pSPARC->dmcomm_phi);

if (pSPARC->Nspden > 1) { // send rho_up, rho_down
D2D(&d2d_sender, &d2d_recvr, gridsizes, pSPARC->DMVertices, pSPARC->electronDens+DMnd, rDMVert,
rho+Nd, pSPARC->dmcomm_phi, sdims, recv_comm, rdims, pSPARC->dmcomm_phi);
D2D(&d2d_sender, &d2d_recvr, gridsizes, pSPARC->DMVertices, pSPARC->electronDens+2*DMnd, rDMVert,
rho+2*Nd, pSPARC->dmcomm_phi, sdims, recv_comm, rdims, pSPARC->dmcomm_phi);
}

// send b
D2D(&d2d_sender, &d2d_recvr, gridsizes, pSPARC->DMVertices, pSPARC->psdChrgDens, rDMVert,
b, pSPARC->dmcomm_phi, sdims, recv_comm, rdims, pSPARC->dmcomm_phi);

// free D2D targets
Free_D2D_Target(&d2d_sender, &d2d_recvr, pSPARC->dmcomm_phi, recv_comm);
if (rank_dmcomm_phi == 0)
MPI_Comm_free(&recv_comm);
} else {
rho = pSPARC->electronDens;
b = pSPARC->psdChrgDens;
}

int n_rho = pSPARC->Nspden/2*2+1;
int length = pSPARC->Nd * n_rho;
if (rank_dmcomm_phi == 0)
{
memcpy(global_rho, rho, length * sizeof(double));
memcpy(global_psdchrgdens, b, pSPARC->Nd * sizeof(double));
}

// free the collected data after printing to file
if (nproc_dmcomm_phi > 1) {
if (rank_dmcomm_phi == 0) {
free(rho);
free(b);
}
}
}

/**
* @brief function to calculate HSMP descriptors
*/
Expand Down
7 changes: 1 addition & 6 deletions src/multipole_features/include/MCSHMainCalc.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@
void Calculate_MCSHDescriptors(SPARC_OBJ *pSPARC, MULTIPOLE_OBJ *mp, const int iterNum);

/**
* @brief function to structure for multipole feature calculation
* @brief function to initialize structure for multipole feature calculation
*/
void Multipole_Initialize(SPARC_OBJ *pSPARC, MULTIPOLE_OBJ *mp);

Expand All @@ -31,8 +31,3 @@ void CalculateHSMPDescriptors(SPARC_OBJ *pSPARC, MULTIPOLE_OBJ *mp, const int it
void CalculateLPMPDescriptors(SPARC_OBJ *pSPARC, MULTIPOLE_OBJ *mp, const int iterNum, double *elecDens, const int commIndex,
const int numParallelComm, const MPI_Comm communicator, int DMVerts[6], const int nFeatures);

/**
* @brief function to collect distributed rho vector into a global vector
*/
void CollectElecDens(SPARC_OBJ *pSPARC, double *global_rho, double *global_psdchrgdens);

0 comments on commit 0242c6d

Please sign in to comment.