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.

diff --git a/src/gromacs/ewald/ewald.cpp b/src/gromacs/ewald/ewald.cpp
index 7b155afde00f20d1a4e0a8e0877f42c6e8ac28e3..5f67078b6805e3bee292a57f3418014e9b68c29e 100644 (file)
--- 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<t_complex, DIM>;
 
-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 7fc31df1e869c996ee9597473d2ad054e500a621..9bf4b684bb43f0a7f3be41b8fcb9c9271eed05d3 100644 (file)
--- a/src/gromacs/ewald/ewald.h
+++ b/src/gromacs/ewald/ewald.h
@@ -68,15 +68,30 @@
 
 #include <stdio.h>
 
+#include <vector>
+
 #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<t_complex> tab_xy;
+    std::vector<t_complex> 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 b329647bf1d5af2f5e90a28e838d1e1f8000ee43..87069cf7ec574a65a1164a36b4e3f557c21def1e 100644 (file)
--- 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 b53ab898c23c567eeaf7a25c6b2172ff45f1f0f1..5bb62cd5b2c378c10c3e407b1d64eae32eba6c8b 100644 (file)
--- 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<int>(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 4ac7bf9141d786f30b2522a048698969a2bad98a..907e79698655124457df87fcf227741571edeabb 100644 (file)
--- 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<gmx_ewald_tab_t>(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 618140e766b2e059fffd4f1b59c5ebbba97fe900..72d5e855df3bfc8121ccae3fc13f57d05e8a5fa8 100644 (file)
--- 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<gmx_ewald_tab_t> ewald_table;
 
     /* Non bonded Parameter lists */
     int               ntype = 0; /* Number of atom types */