From 6fcfe25d8f12fab9a337d645adf9d1da323f2bf0 Mon Sep 17 00:00:00 2001 From: ejjordan Date: Tue, 18 May 2021 20:19:54 +0200 Subject: [PATCH] Remove forcerec and inputrec from ewald_LRcorrection Only pass needed parameters to ewald_LRcorrection, instead of entire forcerec and inputrec. This required also refactoring EwaldBoxZScaler to only take needed parameters instead of a whole inputrec. Also made some names more verbose. --- src/gromacs/ewald/ewald.cpp | 2 +- src/gromacs/ewald/ewald_utils.h | 6 +- src/gromacs/ewald/long_range_correction.cpp | 75 ++++++++------------- src/gromacs/ewald/long_range_correction.h | 15 +++-- src/gromacs/ewald/pme.cpp | 2 +- src/gromacs/ewald/pme_load_balancing.cpp | 2 +- src/gromacs/gmxpreprocess/grompp.cpp | 2 +- src/gromacs/mdlib/force.cpp | 8 ++- 8 files changed, 51 insertions(+), 61 deletions(-) diff --git a/src/gromacs/ewald/ewald.cpp b/src/gromacs/ewald/ewald.cpp index b13e3e7c36..4d04a4f512 100644 --- a/src/gromacs/ewald/ewald.cpp +++ b/src/gromacs/ewald/ewald.cpp @@ -163,7 +163,7 @@ real do_ewald(const t_inputrec& ir, /* Scale box with Ewald wall factor */ matrix scaledBox; - EwaldBoxZScaler boxScaler(ir); + EwaldBoxZScaler boxScaler(inputrecPbcXY2Walls(&ir), ir.wall_ewald_zfac); boxScaler.scaleBox(box, scaledBox); rvec boxDiag; diff --git a/src/gromacs/ewald/ewald_utils.h b/src/gromacs/ewald/ewald_utils.h index 4639d1571b..20db56fceb 100644 --- a/src/gromacs/ewald/ewald_utils.h +++ b/src/gromacs/ewald/ewald_utils.h @@ -100,12 +100,12 @@ public: EwaldBoxZScaler() = delete; /*! \brief Constructor that takes the input record to initialize Ewald box scaling appropriately. */ - EwaldBoxZScaler(const t_inputrec& ir) + EwaldBoxZScaler(bool havePbcXY2Walls, real wallEwaldZfac) { - if (inputrecPbcXY2Walls(&ir)) + if (havePbcXY2Walls) { scaleWithWalls_ = true; - scalingFactor_ = ir.wall_ewald_zfac; + scalingFactor_ = wallEwaldZfac; } else { diff --git a/src/gromacs/ewald/long_range_correction.cpp b/src/gromacs/ewald/long_range_correction.cpp index b5791ead47..b331628599 100644 --- a/src/gromacs/ewald/long_range_correction.cpp +++ b/src/gromacs/ewald/long_range_correction.cpp @@ -65,18 +65,22 @@ * undefined when LJ-PME is not active. This works because * bHaveChargeOrTypePerturbed handles the control flow. */ void ewald_LRcorrection(const int numAtomsLocal, - const t_commrec* cr, + const t_commrec* commrec, int numThreads, int thread, - const t_forcerec& fr, - const t_inputrec& ir, + const real epsilonR, + gmx::ArrayRef qsum, + EwaldGeometry ewaldGeometry, + const real epsilonSurface, + bool havePbcXY2Walls, + real wallEwaldZfac, gmx::ArrayRef chargeA, gmx::ArrayRef chargeB, bool bHaveChargePerturbed, - gmx::ArrayRef x, + gmx::ArrayRef coords, const matrix box, gmx::ArrayRef mu_tot, - gmx::ArrayRef f, + gmx::ArrayRef forces, real* Vcorr_q, real lambda_q, real* dvdlambda_q) @@ -95,10 +99,10 @@ void ewald_LRcorrection(const int numAtomsLocal, /* Scale the Ewald unit cell when dimension z is not periodic */ matrix scaledBox; - EwaldBoxZScaler boxScaler(ir); + EwaldBoxZScaler boxScaler(havePbcXY2Walls, wallEwaldZfac); boxScaler.scaleBox(box, scaledBox); - one_4pi_eps = gmx::c_one4PiEps0 / fr.ic->epsilon_r; + one_4pi_eps = gmx::c_one4PiEps0 / epsilonR; Vexcl_q = 0; dvdl_excl_q = 0; Vdipole[0] = 0; @@ -117,13 +121,13 @@ void ewald_LRcorrection(const int numAtomsLocal, dipole_coeff = 0; real boxVolume = scaledBox[XX][XX] * scaledBox[YY][YY] * scaledBox[ZZ][ZZ]; - switch (ir.ewald_geometry) + switch (ewaldGeometry) { case EwaldGeometry::ThreeD: - if (ir.epsilon_surface != 0) + if (epsilonSurface != 0) { - dipole_coeff = 2 * M_PI * gmx::c_one4PiEps0 - / ((2 * ir.epsilon_surface + fr.ic->epsilon_r) * boxVolume); + dipole_coeff = + 2 * M_PI * gmx::c_one4PiEps0 / ((2 * epsilonSurface + epsilonR) * boxVolume); for (i = 0; (i < DIM); i++) { dipcorrA[i] = 2 * dipole_coeff * mutot[0][i]; @@ -138,19 +142,14 @@ void ewald_LRcorrection(const int numAtomsLocal, for (int q = 0; q < (bHaveChargePerturbed ? 2 : 1); q++) { /* Avoid charge corrections with near-zero net charge */ - if (fabs(fr.qsum[q]) > 1e-4) + if (fabs(qsum[q]) > 1e-4) { - chargecorr[q] = 2 * dipole_coeff * fr.qsum[q]; + chargecorr[q] = 2 * dipole_coeff * qsum[q]; } } break; default: gmx_incons("Unsupported Ewald geometry"); } - if (debug) - { - fprintf(debug, "dipcorr = %8.3f %8.3f %8.3f\n", dipcorrA[XX], dipcorrA[YY], dipcorrA[ZZ]); - fprintf(debug, "mutot = %8.3f %8.3f %8.3f\n", mutot[0][XX], mutot[0][YY], mutot[0][ZZ]); - } const bool bNeedLongRangeCorrection = (dipole_coeff != 0); if (bNeedLongRangeCorrection && !bHaveChargePerturbed) { @@ -158,11 +157,11 @@ void ewald_LRcorrection(const int numAtomsLocal, { for (j = 0; (j < DIM); j++) { - f[i][j] -= dipcorrA[j] * chargeA[i]; + forces[i][j] -= dipcorrA[j] * chargeA[i]; } if (chargecorr[0] != 0) { - f[i][ZZ] += chargecorr[0] * chargeA[i] * x[i][ZZ]; + forces[i][ZZ] += chargecorr[0] * chargeA[i] * coords[i][ZZ]; } } } @@ -172,11 +171,12 @@ void ewald_LRcorrection(const int numAtomsLocal, { for (j = 0; (j < DIM); j++) { - f[i][j] -= L1_q * dipcorrA[j] * chargeA[i] + lambda_q * dipcorrB[j] * chargeB[i]; + forces[i][j] -= L1_q * dipcorrA[j] * chargeA[i] + lambda_q * dipcorrB[j] * chargeB[i]; } if (chargecorr[0] != 0 || chargecorr[1] != 0) { - f[i][ZZ] += (L1_q * chargecorr[0] * chargeA[i] + lambda_q * chargecorr[1]) * x[i][ZZ]; + forces[i][ZZ] += + (L1_q * chargecorr[0] * chargeA[i] + lambda_q * chargecorr[1]) * coords[i][ZZ]; } } } @@ -185,7 +185,7 @@ void ewald_LRcorrection(const int numAtomsLocal, Vself_q[1] = 0; /* Global corrections only on master process */ - if (MASTER(cr) && thread == 0) + if (MASTER(commrec) && thread == 0) { for (q = 0; q < (bHaveChargePerturbed ? 2 : 1); q++) { @@ -195,11 +195,11 @@ void ewald_LRcorrection(const int numAtomsLocal, */ if (dipole_coeff != 0) { - if (ir.ewald_geometry == EwaldGeometry::ThreeD) + if (ewaldGeometry == EwaldGeometry::ThreeD) { Vdipole[q] = dipole_coeff * iprod(mutot[q], mutot[q]); } - else if (ir.ewald_geometry == EwaldGeometry::ThreeDC) + else if (ewaldGeometry == EwaldGeometry::ThreeDC) { Vdipole[q] = dipole_coeff * mutot[q][ZZ] * mutot[q][ZZ]; @@ -215,10 +215,10 @@ void ewald_LRcorrection(const int numAtomsLocal, real sumQZ2 = 0; for (int i = 0; i < numAtomsLocal; i++) { - sumQZ2 += charge[i] * x[i][ZZ] * x[i][ZZ]; + sumQZ2 += charge[i] * coords[i][ZZ] * coords[i][ZZ]; } - Vdipole[q] -= dipole_coeff * fr.qsum[q] - * (sumQZ2 + fr.qsum[q] * box[ZZ][ZZ] * box[ZZ][ZZ] / 12); + Vdipole[q] -= dipole_coeff * qsum[q] + * (sumQZ2 + qsum[q] * box[ZZ][ZZ] * box[ZZ][ZZ] / 12); } } } @@ -233,23 +233,4 @@ void ewald_LRcorrection(const int numAtomsLocal, *Vcorr_q = L1_q * (Vdipole[0] - Vself_q[0]) + lambda_q * (Vdipole[1] - Vself_q[1]) - Vexcl_q; *dvdlambda_q += Vdipole[1] - Vself_q[1] - (Vdipole[0] - Vself_q[0]) - dvdl_excl_q; } - - if (debug) - { - fprintf(debug, "Long Range corrections for Ewald interactions:\n"); - fprintf(debug, - "q2sum = %g, Vself_q=%g\n", - L1_q * fr.q2sum[0] + lambda_q * fr.q2sum[1], - L1_q * Vself_q[0] + lambda_q * Vself_q[1]); - fprintf(debug, "Electrostatic Long Range correction: Vexcl=%g\n", Vexcl_q); - if (MASTER(cr) && thread == 0) - { - if (ir.epsilon_surface > 0 || ir.ewald_geometry == EwaldGeometry::ThreeDC) - { - fprintf(debug, - "Total dipole correction: Vdipole=%g\n", - L1_q * Vdipole[0] + lambda_q * Vdipole[1]); - } - } - } } diff --git a/src/gromacs/ewald/long_range_correction.h b/src/gromacs/ewald/long_range_correction.h index f6fb0c5e8c..b70e3b717d 100644 --- a/src/gromacs/ewald/long_range_correction.h +++ b/src/gromacs/ewald/long_range_correction.h @@ -57,6 +57,7 @@ struct t_commrec; struct t_forcerec; struct t_inputrec; +enum class EwaldGeometry : int; namespace gmx { @@ -69,18 +70,22 @@ class ArrayRef; * Calculate correction for electrostatic surface dipole terms. */ void ewald_LRcorrection(int numAtomsLocal, - const t_commrec* cr, + const t_commrec* commrec, int numThreads, int thread, - const t_forcerec& fr, - const t_inputrec& ir, + real epsilonR, + gmx::ArrayRef qsum, + EwaldGeometry ewaldGeometry, + real epsilonSurface, + bool havePbcXY2Walls, + real wallEwaldZfac, gmx::ArrayRef chargeA, gmx::ArrayRef chargeB, bool bHaveChargePerturbed, - gmx::ArrayRef x, + gmx::ArrayRef coords, const matrix box, gmx::ArrayRef mu_tot, - gmx::ArrayRef f, + gmx::ArrayRef forces, real* Vcorr_q, real lambda_q, real* dvdlambda_q); diff --git a/src/gromacs/ewald/pme.cpp b/src/gromacs/ewald/pme.cpp index 9044c496d2..a3dd835755 100644 --- a/src/gromacs/ewald/pme.cpp +++ b/src/gromacs/ewald/pme.cpp @@ -728,7 +728,7 @@ gmx_pme_t* gmx_pme_init(const t_commrec* cr, // The box requires scaling with nwalls = 2, we store that condition as well // as the scaling factor delete pme->boxScaler; - pme->boxScaler = new EwaldBoxZScaler(*ir); + pme->boxScaler = new EwaldBoxZScaler(inputrecPbcXY2Walls(ir), ir->wall_ewald_zfac); /* If we violate restrictions, generate a fatal error here */ gmx_pme_check_restrictions( diff --git a/src/gromacs/ewald/pme_load_balancing.cpp b/src/gromacs/ewald/pme_load_balancing.cpp index f554e43827..52aaeb9862 100644 --- a/src/gromacs/ewald/pme_load_balancing.cpp +++ b/src/gromacs/ewald/pme_load_balancing.cpp @@ -251,7 +251,7 @@ void pme_loadbal_init(pme_load_balancing_t** pme_lb_p, /* Scale box with Ewald wall factor; note that we pmedata->boxScaler * can't always usedd as it's not available with separate PME ranks. */ - EwaldBoxZScaler boxScaler(ir); + EwaldBoxZScaler boxScaler(inputrecPbcXY2Walls(&ir), ir.wall_ewald_zfac); boxScaler.scaleBox(box, pme_lb->box_start); pme_lb->setup.resize(1); diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index 63f93ae8db..64d6de8c58 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -2343,7 +2343,7 @@ int gmx_grompp(int argc, char* argv[]) { /* Calculate the optimal grid dimensions */ matrix scaledBox; - EwaldBoxZScaler boxScaler(*ir); + EwaldBoxZScaler boxScaler(inputrecPbcXY2Walls(ir), ir->wall_ewald_zfac); boxScaler.scaleBox(state.box, scaledBox); if (ir->nkx > 0 && ir->nky > 0 && ir->nkz > 0) diff --git a/src/gromacs/mdlib/force.cpp b/src/gromacs/mdlib/force.cpp index 8f07ca6231..1112344817 100644 --- a/src/gromacs/mdlib/force.cpp +++ b/src/gromacs/mdlib/force.cpp @@ -165,8 +165,12 @@ void calculateLongRangeNonbondeds(t_forcerec* fr, cr, nthreads, t, - *fr, - ir, + fr->ic->epsilon_r, + fr->qsum, + ir.ewald_geometry, + ir.epsilon_surface, + inputrecPbcXY2Walls(&ir), + ir.wall_ewald_zfac, md->chargeA ? gmx::constArrayRefFromArray(md->chargeA, md->nr) : gmx::ArrayRef{}, md->chargeB ? gmx::constArrayRefFromArray(md->chargeB, md->nr) -- 2.22.0