From 8e655583a6fd9bd2c027f4b92b9921dba5ea9fe6 Mon Sep 17 00:00:00 2001 From: Joe Jordan Date: Mon, 22 Feb 2021 19:00:38 +0000 Subject: [PATCH] Hold unique_ptr to gmx_ewald_tab_t in forcerec Proper constructor and destructor are added. The field eir, which was anyway recomputed with every call to do_ewald, is removed from the struct and made local to do_ewald. Since the type cvec is only used by eir, it is moved to reside in ewald.cpp. --- src/gromacs/ewald/ewald.cpp | 49 ++++++++++++++-------------------- src/gromacs/ewald/ewald.h | 19 +++++++++++-- src/gromacs/math/gmxcomplex.h | 4 +-- src/gromacs/mdlib/force.cpp | 2 +- src/gromacs/mdlib/forcerec.cpp | 2 +- src/gromacs/mdtypes/forcerec.h | 2 +- 6 files changed, 41 insertions(+), 37 deletions(-) diff --git a/src/gromacs/ewald/ewald.cpp b/src/gromacs/ewald/ewald.cpp index 7b155afde0..5f67078b68 100644 --- a/src/gromacs/ewald/ewald.cpp +++ b/src/gromacs/ewald/ewald.cpp @@ -71,30 +71,23 @@ #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/smalloc.h" -struct gmx_ewald_tab_t -{ - int nx, ny, nz, kmax; - cvec** eir; - t_complex *tab_xy, *tab_qxyz; -}; +using cvec = std::array; -void init_ewald_tab(struct gmx_ewald_tab_t** et, const t_inputrec& ir, FILE* fp) +gmx_ewald_tab_t::gmx_ewald_tab_t(const t_inputrec& ir, FILE* fp) { - snew(*et, 1); if (fp) { fprintf(fp, "Will do ordinary reciprocal space Ewald sum.\n"); } - (*et)->nx = ir.nkx + 1; - (*et)->ny = ir.nky + 1; - (*et)->nz = ir.nkz + 1; - (*et)->kmax = std::max((*et)->nx, std::max((*et)->ny, (*et)->nz)); - (*et)->eir = nullptr; - (*et)->tab_xy = nullptr; - (*et)->tab_qxyz = nullptr; + nx = ir.nkx + 1; + ny = ir.nky + 1; + nz = ir.nkz + 1; + kmax = std::max(nx, std::max(ny, nz)); } +gmx_ewald_tab_t::~gmx_ewald_tab_t() = default; + //! Calculates wave vectors. static void calc_lll(const rvec box, rvec lll) { @@ -158,6 +151,7 @@ real do_ewald(const t_inputrec& ir, int lowiy, lowiz, ix, iy, iz, n, q; real tmp, cs, ss, ak, akv, mx, my, mz, m2, scale; gmx_bool bFreeEnergy; + cvec** eir; if (cr != nullptr) { @@ -181,23 +175,20 @@ real do_ewald(const t_inputrec& ir, /* 1/(Vol*e0) */ real scaleRecip = 4.0 * M_PI / (boxDiag[XX] * boxDiag[YY] * boxDiag[ZZ]) * ONE_4PI_EPS0 / ir.epsilon_r; - if (!et->eir) /* allocate if we need to */ + snew(eir, et->kmax); + for (n = 0; n < et->kmax; n++) { - snew(et->eir, et->kmax); - for (n = 0; n < et->kmax; n++) - { - snew(et->eir[n], natoms); - } - snew(et->tab_xy, natoms); - snew(et->tab_qxyz, natoms); + snew(eir[n], natoms); } + et->tab_xy.resize(natoms); + et->tab_qxyz.resize(natoms); bFreeEnergy = (ir.efep != FreeEnergyPerturbationType::No); clear_mat(lrvir); calc_lll(boxDiag, lll); - tabulateStructureFactors(natoms, x, et->kmax, et->eir, lll); + tabulateStructureFactors(natoms, x, et->kmax, eir, lll); for (q = 0; q < (bFreeEnergy ? 2 : 1); q++) { @@ -229,14 +220,14 @@ real do_ewald(const t_inputrec& ir, { for (n = 0; n < natoms; n++) { - et->tab_xy[n] = cmul(et->eir[ix][n][XX], et->eir[iy][n][YY]); + et->tab_xy[n] = cmul(eir[ix][n][XX], eir[iy][n][YY]); } } else { for (n = 0; n < natoms; n++) { - et->tab_xy[n] = cmul(et->eir[ix][n][XX], conjugate(et->eir[-iy][n][YY])); + et->tab_xy[n] = cmul(eir[ix][n][XX], conjugate(eir[-iy][n][YY])); } } for (iz = lowiz; iz < et->nz; iz++) @@ -249,15 +240,15 @@ real do_ewald(const t_inputrec& ir, { for (n = 0; n < natoms; n++) { - et->tab_qxyz[n] = rcmul(charge[n], cmul(et->tab_xy[n], et->eir[iz][n][ZZ])); + et->tab_qxyz[n] = rcmul(charge[n], cmul(et->tab_xy[n], eir[iz][n][ZZ])); } } else { for (n = 0; n < natoms; n++) { - et->tab_qxyz[n] = rcmul( - charge[n], cmul(et->tab_xy[n], conjugate(et->eir[-iz][n][ZZ]))); + et->tab_qxyz[n] = + rcmul(charge[n], cmul(et->tab_xy[n], conjugate(eir[-iz][n][ZZ]))); } } diff --git a/src/gromacs/ewald/ewald.h b/src/gromacs/ewald/ewald.h index 7fc31df1e8..9bf4b684bb 100644 --- a/src/gromacs/ewald/ewald.h +++ b/src/gromacs/ewald/ewald.h @@ -68,15 +68,30 @@ #include +#include + #include "gromacs/math/vectypes.h" #include "gromacs/utility/real.h" struct t_commrec; struct t_forcerec; struct t_inputrec; +struct t_complex; + +struct gmx_ewald_tab_t +{ + gmx_ewald_tab_t(const t_inputrec& ir, FILE* fp); + + ~gmx_ewald_tab_t(); + + int nx; + int ny; + int nz; + int kmax; -/* Forward declaration of type for managing Ewald tables */ -struct gmx_ewald_tab_t; + std::vector tab_xy; + std::vector tab_qxyz; +}; /*! \brief Initialize the tables used in the Ewald long-ranged part */ void init_ewald_tab(struct gmx_ewald_tab_t** et, const t_inputrec& ir, FILE* fp); diff --git a/src/gromacs/math/gmxcomplex.h b/src/gromacs/math/gmxcomplex.h index b329647bf1..87069cf7ec 100644 --- a/src/gromacs/math/gmxcomplex.h +++ b/src/gromacs/math/gmxcomplex.h @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2012,2014,2017,2018,2019, by the GROMACS development team, led by + * Copyright (c) 2012,2014,2017,2018,2019,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -47,8 +47,6 @@ struct t_complex real re, im; }; -typedef t_complex cvec[DIM]; - static t_complex rcmul(real r, t_complex c) { t_complex d; diff --git a/src/gromacs/mdlib/force.cpp b/src/gromacs/mdlib/force.cpp index b53ab898c2..5bb62cd5b2 100644 --- a/src/gromacs/mdlib/force.cpp +++ b/src/gromacs/mdlib/force.cpp @@ -282,7 +282,7 @@ void calculateLongRangeNonbondeds(t_forcerec* fr, fr->ic->ewaldcoeff_q, lambda[static_cast(FreeEnergyPerturbationCouplingType::Coul)], &ewaldOutput.dvdl[FreeEnergyPerturbationCouplingType::Coul], - fr->ewald_table); + fr->ewald_table.get()); } /* Note that with separate PME nodes we get the real energies later */ diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp index 4ac7bf9141..907e796986 100644 --- a/src/gromacs/mdlib/forcerec.cpp +++ b/src/gromacs/mdlib/forcerec.cpp @@ -1107,7 +1107,7 @@ void init_forcerec(FILE* fp, /* TODO: Replace this Ewald table or move it into interaction_const_t */ if (ir.coulombtype == CoulombInteractionType::Ewald) { - init_ewald_tab(&(fr->ewald_table), ir, fp); + fr->ewald_table = std::make_unique(ir, fp); } /* Electrostatics: Translate from interaction-setting-in-mdp-file to kernel interaction format */ diff --git a/src/gromacs/mdtypes/forcerec.h b/src/gromacs/mdtypes/forcerec.h index 618140e766..72d5e855df 100644 --- a/src/gromacs/mdtypes/forcerec.h +++ b/src/gromacs/mdtypes/forcerec.h @@ -265,7 +265,7 @@ struct t_forcerec LongRangeVdW ljpme_combination_rule = LongRangeVdW::Geom; /* PME/Ewald stuff */ - struct gmx_ewald_tab_t* ewald_table = nullptr; + std::unique_ptr ewald_table; /* Non bonded Parameter lists */ int ntype = 0; /* Number of atom types */ -- 2.22.0