From 037839d0e598f12ee424520ca5507e1fa8e604a8 Mon Sep 17 00:00:00 2001 From: ejjordan Date: Wed, 19 May 2021 18:28:07 +0200 Subject: [PATCH] Remove inputrec from do_ewald Also rename some parameters more verbosely and remove a dead code block. --- src/gromacs/ewald/ewald.cpp | 35 ++++++++++++++++------------------- src/gromacs/ewald/ewald.h | 12 ++++++++---- src/gromacs/mdlib/force.cpp | 5 ++++- 3 files changed, 28 insertions(+), 24 deletions(-) diff --git a/src/gromacs/ewald/ewald.cpp b/src/gromacs/ewald/ewald.cpp index 4d04a4f512..65e31644c3 100644 --- a/src/gromacs/ewald/ewald.cpp +++ b/src/gromacs/ewald/ewald.cpp @@ -131,13 +131,16 @@ static void tabulateStructureFactors(int natom, gmx::ArrayRef x } } -real do_ewald(const t_inputrec& ir, - gmx::ArrayRef x, - gmx::ArrayRef f, +real do_ewald(bool havePbcXY2Walls, + real wallEwaldZfac, + real epsilonR, + FreeEnergyPerturbationType freeEnergyPerturbationType, + gmx::ArrayRef coords, + gmx::ArrayRef forces, gmx::ArrayRef chargeA, gmx::ArrayRef chargeB, const matrix box, - const t_commrec* cr, + const t_commrec* commrec, int natoms, matrix lrvir, real ewaldcoeff, @@ -153,9 +156,9 @@ real do_ewald(const t_inputrec& ir, cvec** eir; bool bFreeEnergy; - if (cr != nullptr) + if (commrec != nullptr) { - if (PAR(cr)) + if (PAR(commrec)) { gmx_fatal(FARGS, "No parallel Ewald. Use PME instead.\n"); } @@ -163,7 +166,7 @@ real do_ewald(const t_inputrec& ir, /* Scale box with Ewald wall factor */ matrix scaledBox; - EwaldBoxZScaler boxScaler(inputrecPbcXY2Walls(&ir), ir.wall_ewald_zfac); + EwaldBoxZScaler boxScaler(havePbcXY2Walls, wallEwaldZfac); boxScaler.scaleBox(box, scaledBox); rvec boxDiag; @@ -173,8 +176,7 @@ real do_ewald(const t_inputrec& ir, } /* 1/(Vol*e0) */ - real scaleRecip = - 4.0 * M_PI / (boxDiag[XX] * boxDiag[YY] * boxDiag[ZZ]) * gmx::c_one4PiEps0 / ir.epsilon_r; + real scaleRecip = 4.0 * M_PI / (boxDiag[XX] * boxDiag[YY] * boxDiag[ZZ]) * gmx::c_one4PiEps0 / epsilonR; snew(eir, et->kmax); for (n = 0; n < et->kmax; n++) @@ -184,12 +186,12 @@ real do_ewald(const t_inputrec& ir, et->tab_xy.resize(natoms); et->tab_qxyz.resize(natoms); - bFreeEnergy = (ir.efep != FreeEnergyPerturbationType::No); + bFreeEnergy = (freeEnergyPerturbationType != FreeEnergyPerturbationType::No); clear_mat(lrvir); calc_lll(boxDiag, lll); - tabulateStructureFactors(natoms, x, et->kmax, eir, lll); + tabulateStructureFactors(natoms, coords, et->kmax, eir, lll); gmx::ArrayRef charge; for (q = 0; q < (bFreeEnergy ? 2 : 1); q++) @@ -272,14 +274,9 @@ real do_ewald(const t_inputrec& ir, { /*tmp=scale*ak*(cs*tab_qxyz[n].im-ss*tab_qxyz[n].re);*/ tmp = scale * ak * (cs * et->tab_qxyz[n].im - ss * et->tab_qxyz[n].re); - f[n][XX] += tmp * mx * 2 * scaleRecip; - f[n][YY] += tmp * my * 2 * scaleRecip; - f[n][ZZ] += tmp * mz * 2 * scaleRecip; -#if 0 - f[n][XX] += tmp*mx; - f[n][YY] += tmp*my; - f[n][ZZ] += tmp*mz; -#endif + forces[n][XX] += tmp * mx * 2 * scaleRecip; + forces[n][YY] += tmp * my * 2 * scaleRecip; + forces[n][ZZ] += tmp * mz * 2 * scaleRecip; } lowiz = 1 - et->nz; } diff --git a/src/gromacs/ewald/ewald.h b/src/gromacs/ewald/ewald.h index cfb6cd2a86..809bc8f4f7 100644 --- a/src/gromacs/ewald/ewald.h +++ b/src/gromacs/ewald/ewald.h @@ -77,6 +77,7 @@ struct t_commrec; struct t_forcerec; struct t_inputrec; struct t_complex; +enum class FreeEnergyPerturbationType : int; namespace gmx { @@ -100,13 +101,16 @@ struct gmx_ewald_tab_t }; /*! \brief Do the long-ranged part of an Ewald calculation */ -real do_ewald(const t_inputrec& ir, - gmx::ArrayRef x, - gmx::ArrayRef f, +real do_ewald(bool havePbcXY2Walls, + real wallEwaldZfac, + real epsilonR, + FreeEnergyPerturbationType freeEnergyPerturbationType, + gmx::ArrayRef coords, + gmx::ArrayRef forces, gmx::ArrayRef chargeA, gmx::ArrayRef chargeB, const matrix box, - const t_commrec* cr, + const t_commrec* commrec, int natoms, matrix lrvir, real ewaldcoeff, diff --git a/src/gromacs/mdlib/force.cpp b/src/gromacs/mdlib/force.cpp index 1112344817..3ac2484a2b 100644 --- a/src/gromacs/mdlib/force.cpp +++ b/src/gromacs/mdlib/force.cpp @@ -281,7 +281,10 @@ void calculateLongRangeNonbondeds(t_forcerec* fr, if (fr->ic->eeltype == CoulombInteractionType::Ewald) { - Vlr_q = do_ewald(ir, + Vlr_q = do_ewald(inputrecPbcXY2Walls(&ir), + ir.wall_ewald_zfac, + ir.epsilon_r, + ir.efep, coordinates, forceWithVirial->force_, gmx::arrayRefFromArray(md->chargeA, md->nr), -- 2.22.0