Refactor reference PME (#5131)

platforms/cpu/src/CpuNonbondedForce.cpp

@@ -255,28 +255,26 @@ void CpuNonbondedForce::calculateReciprocalIxn(int numberOfAtoms, float* posq, c
     float recipCoeff               = (float)(ONE_4PI_EPS0*4*PI_M/(periodicBoxVectors[0][0] * periodicBoxVectors[1][1] * periodicBoxVectors[2][2]) /epsilon);
 
     if (pme) {
-        pme_t pmedata;
-        pme_init(&pmedata, alphaEwald, numberOfAtoms, meshDim, 5, 1);
+        ReferencePME pme(alphaEwald, numberOfAtoms, meshDim, 5, 1);
         vector<double> charges(numberOfAtoms);
         for (int i = 0; i < numberOfAtoms; i++)
             charges[i] = posq[4*i+3];
         double recipEnergy = 0.0;
-        pme_exec(pmedata, atomCoordinates, forces, charges, periodicBoxVectors, &recipEnergy);
+        pme.exec(atomCoordinates, forces, charges, periodicBoxVectors, recipEnergy);
         if (totalEnergy)
             *totalEnergy += recipEnergy;
-        pme_destroy(pmedata);
 
         if (ljpme) {
             // Dispersion reciprocal space terms
-            pme_init(&pmedata,alphaDispersionEwald,numberOfAtoms,dispersionMeshDim,5,1);
+            ReferencePME ljpme(alphaDispersionEwald,numberOfAtoms,dispersionMeshDim,5,1);
 
-            std::vector<Vec3> dpmeforces;
+            vector<Vec3> dpmeforces;
             for (int i = 0; i < numberOfAtoms; i++){
                 charges[i] = C6params[i];
                 dpmeforces.push_back(Vec3());
             }
             double recipDispersionEnergy = 0.0;
-            pme_exec_dpme(pmedata,atomCoordinates,dpmeforces,charges,periodicBoxVectors,&recipDispersionEnergy);
+            ljpme.exec_dpme(atomCoordinates, dpmeforces, charges, periodicBoxVectors, recipDispersionEnergy);
             for (int i = 0; i < numberOfAtoms; i++){
                 forces[i][0] += dpmeforces[i][0];
                 forces[i][1] += dpmeforces[i][1];
@@ -284,8 +282,6 @@ void CpuNonbondedForce::calculateReciprocalIxn(int numberOfAtoms, float* posq, c
             }
             if (totalEnergy)
                 *totalEnergy += recipDispersionEnergy;
-
-            pme_destroy(pmedata);
         }
 
     }

platforms/reference/include/ReferenceConstantPotential.h

@@ -90,14 +90,14 @@ public:
      * @param elecToSys              mapping from electrode particle indices to system particle indices
      * @param electrodeParamArray    electrode particle parameters
      * @param conp                   constant potential derivative evaluation class
-     * @param pmeData                reference PME solver
+     * @param pme                    reference PME solver
      */
     virtual void solve(int numParticles, int numElectrodeParticles,
         const std::vector<Vec3>& posData, std::vector<double>& charges,
         const std::vector<std::set<int>>& exclusions,
         const std::vector<int>& sysToElec, const std::vector<int>& elecToSys,
         const std::vector<std::array<double, 3> >& electrodeParamArray,
-        ReferenceConstantPotential& conp, pme_t pmeData) = 0;
+        ReferenceConstantPotential& conp, ReferencePME& pme) = 0;
 };
 
 /**
@@ -149,14 +149,14 @@ public:
      * @param elecToSys              mapping from electrode particle indices to system particle indices
      * @param electrodeParamArray    electrode particle parameters
      * @param conp                   constant potential derivative evaluation class
-     * @param pmeData                reference PME solver
+     * @param pme                    reference PME solver
      */
     void solve(int numParticles, int numElectrodeParticles,
         const std::vector<Vec3>& posData, std::vector<double>& charges,
         const std::vector<std::set<int>>& exclusions,
         const std::vector<int>& sysToElec, const std::vector<int>& elecToSys,
         const std::vector<std::array<double, 3> >& electrodeParamArray,
-        ReferenceConstantPotential& conp, pme_t pmeData);
+        ReferenceConstantPotential& conp, ReferencePME& pme);
 };
 
 /**
@@ -208,14 +208,14 @@ public:
      * @param elecToSys              mapping from electrode particle indices to system particle indices
      * @param electrodeParamArray    electrode particle parameters
      * @param conp                   constant potential derivative evaluation class
-     * @param pmeData                reference PME solver
+     * @param pme                    reference PME solver
      */
     void solve(int numParticles, int numElectrodeParticles,
         const std::vector<Vec3>& posData, std::vector<double>& charges,
         const std::vector<std::set<int>>& exclusions,
         const std::vector<int>& sysToElec, const std::vector<int>& elecToSys,
         const std::vector<std::array<double, 3> >& electrodeParamArray,
-        ReferenceConstantPotential& conp, pme_t pmeData);
+        ReferenceConstantPotential& conp, ReferencePME& pme);
 };
 
 /**
@@ -313,7 +313,7 @@ private:
      * @param elecToSys              mapping from electrode particle indices to system particle indices
      * @param electrodeParamArray    electrode particle parameters
      * @param energy                 output system energy
-     * @param pmeData                reference PME solver
+     * @param pme                    reference PME solver
      */
     void getEnergyForces(int numParticles, int numElectrodeParticles,
         const std::vector<Vec3>& posData, std::vector<Vec3>& forceData,
@@ -321,7 +321,7 @@ private:
         const std::vector<std::set<int>>& exclusions,
         const std::vector<int>& sysToElec, const std::vector<int>& elecToSys,
         const std::vector<std::array<double, 3> >& electrodeParamArray,
-        double* energy, pme_t pmeData);
+        double* energy, ReferencePME& pme);
     /**
      * Computes energy derivatives with respect to charges.
      * 
@@ -334,14 +334,14 @@ private:
      * @param elecToSys              mapping from electrode particle indices to system particle indices
      * @param electrodeParamArray    electrode particle parameters
      * @param chargeDerivatives      output charge derivatives
-     * @param pmeData                reference PME solver
+     * @param pme                    reference PME solver
      */
     void getDerivatives(int numParticles, int numElectrodeParticles,
         const std::vector<Vec3>& posData, std::vector<double>& charges,
         const std::vector<std::set<int>>& exclusions,
         const std::vector<int>& sysToElec, const std::vector<int>& elecToSys,
         const std::vector<std::array<double, 3> >& electrodeParamArray,
-        std::vector<double>& chargeDerivatives, pme_t pmeData);
+        std::vector<double>& chargeDerivatives, ReferencePME& pme);
 };
 
 } // namespace OpenMM

platforms/reference/include/ReferencePME.h

 /* 
  * Reference implementation of PME reciprocal space interactions.
  * 
- * Copyright (c) 2009, Erik Lindahl, Rossen Apostolov, Szilard Pall
+ * Copyright (c) 2009-2025 Erik Lindahl, Rossen Apostolov, Szilard Pall, Peter Eastman
  * All rights reserved.
  * Contact: lindahl@cbr.su.se Stockholm University, Sweden.
  * 
@@ -34,22 +34,19 @@
 
 #include "openmm/Vec3.h"
 #include "openmm/internal/windowsExport.h"
+#include <array>
+#include <complex>
 #include <vector>
 
 namespace OpenMM {
 
-typedef double rvec[3];
-
-
-typedef struct pme *
-pme_t;
-
-/*
+class OPENMM_EXPORT ReferencePME {
+public:
+    /*
      * Initialize a PME calculation and set up data structures
      *
      * Arguments:
      * 
- * ppme        Pointer to an opaque pme_t object
      * ewaldcoeff  Coefficient derived from the beta factor to participate
      *             direct/reciprocal space. See gromacs code for documentation!
      *             We assume that you are using nm units...
@@ -58,77 +55,96 @@ pme_t;
      * pme_order   Interpolation order, almost always 4
      * epsilon_r   Dielectric coefficient, typically 1.0.
      */
-int OPENMM_EXPORT
-pme_init(pme_t* ppme,
-         double ewaldcoeff,
-         int natoms,
-         const int ngrid[3],
-         int pme_order,
-         double epsilon_r);
+    ReferencePME(double ewaldcoeff, int natoms, const int ngrid[3], int pme_order, double epsilon_r);
 
-/*
+    /*
      * Evaluate reciprocal space PME energy and forces.
      *
      * Args:
      *
- * pme                 Opaque pme_t object, must have been initialized with pme_init()
      * atomCoordinates     Pointer to coordinate data array (nm)
      * forces              Pointer to force data array (will be written as kJ/mol/nm)
      * charges             Array of charges (units of e)
      * periodicBoxVectors  Simulation cell dimensions (nm)
      * energy              Total energy (will be written in units of kJ/mol)
      */
-int OPENMM_EXPORT
-pme_exec(pme_t pme,
-         const std::vector<OpenMM::Vec3>& atomCoordinates,
-         std::vector<OpenMM::Vec3>& forces,
-         const std::vector<double>& charges,
-         const OpenMM::Vec3 periodicBoxVectors[3],
-         double* energy);
+    void exec(const std::vector<OpenMM::Vec3>& atomCoordinates, std::vector<OpenMM::Vec3>& forces, const std::vector<double>& charges,
+            const OpenMM::Vec3 periodicBoxVectors[3], double& energy);
 
-/*
+    /*
      * Evaluate reciprocal space PME energy and charge derivatives.
      *
      * Args:
      *
- * pme                 Opaque pme_t object, must have been initialized with pme_init()
      * atomCoordinates     Pointer to coordinate data array (nm)
      * chargeDerivatives   Pointer to charge derivative data array (will be written as kJ/mol/e)
      * chargeIndices       Pointer to array of indices of particles to compute charge derivatives for
      * charges             Array of charges (units of e)
      * periodicBoxVectors  Simulation cell dimensions (nm)
      */
-int OPENMM_EXPORT
-pme_exec_charge_derivatives(pme_t pme,
-                            const std::vector<OpenMM::Vec3>& atomCoordinates,
-                            std::vector<double>& chargeDerivatives,
-                            const std::vector<int>& chargeIndices,
-                            const std::vector<double>& charges,
-                            const OpenMM::Vec3 periodicBoxVectors[3]);
+    void exec_charge_derivatives(const std::vector<OpenMM::Vec3>& atomCoordinates, std::vector<double>& chargeDerivatives,
+            const std::vector<int>& chargeIndices, const std::vector<double>& charges, const OpenMM::Vec3 periodicBoxVectors[3]);
 
-/**
+    /**
      * Evaluate reciprocal space PME dispersion energy and forces.
      *
      * Args:
      *
- * pme                 Opaque pme_t object, must have been initialized with pme_init()
      * atomCoordinates     Pointer to coordinate data array (nm)
      * forces              Pointer to force data array (will be written as kJ/mol/nm)
      * c6s                 Array of c6 coefficients (units of sqrt(kJ/mol).nm^3 )
      * periodicBoxVectors  Simulation cell dimensions (nm)
      * energy              Total energy (will be written in units of kJ/mol)
      */
-int OPENMM_EXPORT
-pme_exec_dpme(pme_t pme,
-              const std::vector<OpenMM::Vec3>& atomCoordinates,
-              std::vector<OpenMM::Vec3>& forces,
-              const std::vector<double>& c6s,
-              const OpenMM::Vec3 periodicBoxVectors[3],
-              double* energy);
+    void exec_dpme(const std::vector<OpenMM::Vec3>& atomCoordinates, std::vector<OpenMM::Vec3>& forces, const std::vector<double>& c6s,
+            const OpenMM::Vec3 periodicBoxVectors[3], double& energy);
+private:
+    void calculate_bsplines_moduli();
+    void update_grid_index_and_fraction(const std::vector<Vec3>& atomCoordinates, const Vec3 recipBoxVectors[3]);
+    void update_bsplines();
+    void grid_spread_charge(const std::vector<double>& charges);
+    void pme_reciprocal_convolution(const Vec3 periodicBoxVectors[3], const Vec3 recipBoxVectors[3], double& energy);
+    void dpme_reciprocal_convolution(const Vec3 periodicBoxVectors[3], const Vec3 recipBoxVectors[3], double& energy);
+    void grid_interpolate_force(const Vec3 recipBoxVectors[3], const std::vector<double>& charges, std::vector<Vec3>& forces);
+    void grid_interpolate_charge_derivatives(const Vec3 recipBoxVectors[3], const std::vector<double>& charges,
+            std::vector<double>& chargeDerivatives, const std::vector<int>& chargeIndices);
+    int natoms;
+    double ewaldcoeff;
+    std::vector<std::complex<double> > grid; /* Memory for the grid we spread charges on.
+                                              * Element (i,j,k) is accessed as:
+                                              * grid[i*ngrid[1]*ngrid[2] + j*ngrid[2] + k] */
+    int ngrid[3];             /* Total grid dimensions (all data is complex!) */
+    int order;                /* PME interpolation order. Almost always 4 */
+
+    /* Data for bspline interpolation, see the Essman PME paper */
+    std::vector<double> bsplines_moduli[3];   /* 3 pointers, to x/y/z bspline moduli, each of length ngrid[x/y/z]   */
+    std::vector<double> bsplines_theta[3];    /* each of x/y/z has length order*natoms */
+    std::vector<double> bsplines_dtheta[3];   /* each of x/y/z has length order*natoms */
+
+    std::vector<std::array<int, 3> > particleindex; /* Array of length natoms. Each element is
+                                                     * three ints that specify the grid
+                                                     * indices for that particular atom. Updated every step! */
+    std::vector<Vec3> particlefraction; /* Array of length natoms. Fractional offset in the grid for
+                                         * each atom in all three dimensions. */
+
+    /* Further explanation of index/fraction:
+     *
+     * Assume we have a cell of size 10*10*10nm, and a total grid dimension of 100*100*100 cells.
+     * In other words, each cell is 0.1*0.1*0.1 nm.
+     *
+     * If particle i has coordinates { 0.543 , 6.235 , -0.73 }, we will get:
+     *
+     * particleindex[i]    = { 5 , 62 , 92 }         (-0.73 + 10 = 9.27, we always apply PBC for grid calculations!)
+     * particlefraction[i] = { 0.43 , 0.35 , 0.7 }   (this is the fraction of the cell length where the atom is)
+     *
+     * (The reason for precaculating / storing these is that it gets a bit more complex for triclinic cells :-)
+     *
+     * In the current code version we might assume that a particle is not more than a whole box length away from
+     * the central cell, i.e., in this case we would assume all coordinates fall in -10 nm < x,y,z < 20 nm.
+     */
 
-/* Release all memory in pme structure */
-int OPENMM_EXPORT
-pme_destroy(pme_t    pme);
+    double epsilon_r; /* Dielectric coefficient to use, typically 1.0 */
+};
 
 } // namespace OpenMM
 

platforms/reference/src/SimTKReference/ReferenceConstantPotential.cpp

@@ -29,7 +29,6 @@
 
 #include "ReferenceConstantPotential.h"
 #include "ReferenceForce.h"
-#include "ReferencePME.h"
 #include "SimTKOpenMMUtilities.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/MSVC_erfc.h"
@@ -97,22 +96,21 @@ void ReferenceConstantPotentialMatrixSolver::update(
     std::vector<double> dUdQ0(numElectrodeParticles);
     std::vector<double> dUdQ(numElectrodeParticles);
     
-    pme_t pmeData;
-    pme_init(&pmeData, conp.ewaldAlpha, numParticles, conp.gridSize, 5, 1);
+    ReferencePME pme(conp.ewaldAlpha, numParticles, conp.gridSize, 5, 1);
 
     // Get derivatives when all electrode charges are zeroed.
     std::vector<double> q(charges);
     for (int ii = 0; ii < numElectrodeParticles; ii++) {
         q[elecToSys[ii]] = 0.0;
     }
-    conp.getDerivatives(numParticles, numElectrodeParticles, posData, q, exclusions, sysToElec, elecToSys, electrodeParamArray, dUdQ0, pmeData);
+    conp.getDerivatives(numParticles, numElectrodeParticles, posData, q, exclusions, sysToElec, elecToSys, electrodeParamArray, dUdQ0, pme);
 
     for (int ii = 0; ii < numElectrodeParticles; ii++) {
         int i = elecToSys[ii];
 
         // Get derivatives when one electrode charge is set.
         q[i] = 1.0;
-        conp.getDerivatives(numParticles, numElectrodeParticles, posData, q, exclusions, sysToElec, elecToSys, electrodeParamArray, dUdQ, pmeData);
+        conp.getDerivatives(numParticles, numElectrodeParticles, posData, q, exclusions, sysToElec, elecToSys, electrodeParamArray, dUdQ, pme);
         q[i] = 0.0;
 
         // Set matrix elements, subtracting zero charge derivatives so that the
@@ -123,8 +121,6 @@ void ReferenceConstantPotentialMatrixSolver::update(
         A[ii][ii] = dUdQ[ii] - dUdQ0[ii];
     }
 
-    pme_destroy(pmeData);
-
     // Compute Cholesky decomposition representation of the inverse.
     capacitance = JAMA::Cholesky<double>(A);
     if (!capacitance.is_spd()) {
@@ -162,7 +158,7 @@ void ReferenceConstantPotentialMatrixSolver::solve(
     const std::vector<int>& elecToSys,
     const std::vector<std::array<double, 3> >& electrodeParamArray,
     ReferenceConstantPotential& conp,
-    pme_t pmeData
+    ReferencePME& pme
 ) {
     // Solves for charges using the matrix method.
 
@@ -171,7 +167,7 @@ void ReferenceConstantPotentialMatrixSolver::solve(
         charges[elecToSys[ii]] = 0.0;
     }
     std::vector<double> b(numElectrodeParticles);
-    conp.getDerivatives(numParticles, numElectrodeParticles, posData, charges, exclusions, sysToElec, elecToSys, electrodeParamArray, b, pmeData);
+    conp.getDerivatives(numParticles, numElectrodeParticles, posData, charges, exclusions, sysToElec, elecToSys, electrodeParamArray, b, pme);
     for (int ii = 0; ii < numElectrodeParticles; ii++) {
         b[ii] = -b[ii];
     }
@@ -236,14 +232,12 @@ void ReferenceConstantPotentialCGSolver::update(
         // calculation.  This will actually vary slightly with position but only due
         // to finite accuracy of the PME splines, so it is fine to assume it will be
         // constant for the preconditioner.
-        pme_t pmeData;
-        pme_init(&pmeData, conp.ewaldAlpha, 1, conp.gridSize, 5, 1);
+        ReferencePME pme(conp.ewaldAlpha, 1, conp.gridSize, 5, 1);
         std::vector<Vec3> pmePosData(1);
         std::vector<double> pmeChargeDerivatives(1);
         std::vector<int> pmeElectrodeIndices(1);
         std::vector<double> pmeCharges(1, 1.0);
-        pme_exec_charge_derivatives(pmeData, pmePosData, pmeChargeDerivatives, pmeElectrodeIndices, pmeCharges, boxVectors);
-        pme_destroy(pmeData);
+        pme.exec_charge_derivatives(pmePosData, pmeChargeDerivatives, pmeElectrodeIndices, pmeCharges, boxVectors);
 
         // The diagonal has a contribution from reciprocal space, Ewald
         // self-interaction, Ewald neutralizing plasma, Gaussian self-interaction,
@@ -274,7 +268,7 @@ void ReferenceConstantPotentialCGSolver::solve(
     const std::vector<int>& elecToSys,
     const std::vector<std::array<double, 3> >& electrodeParamArray,
     ReferenceConstantPotential& conp,
-    pme_t pmeData
+    ReferencePME& pme
 ) {
     // Solves for charges using the conjugate gradient method.
     
@@ -302,7 +296,7 @@ void ReferenceConstantPotentialCGSolver::solve(
     }
 
     // Evaluate the initial gradient Aq - b.
-    conp.getDerivatives(numParticles, numElectrodeParticles, posData, charges, exclusions, sysToElec, elecToSys, electrodeParamArray, grad, pmeData);
+    conp.getDerivatives(numParticles, numElectrodeParticles, posData, charges, exclusions, sysToElec, elecToSys, electrodeParamArray, grad, pme);
 
     // Project the initial gradient without preconditioning.
     offset = 0.0;
@@ -332,7 +326,7 @@ void ReferenceConstantPotentialCGSolver::solve(
         q[ii] = charges[i];
         charges[i] = 0.0;
     }
-    conp.getDerivatives(numParticles, numElectrodeParticles, posData, charges, exclusions, sysToElec, elecToSys, electrodeParamArray, grad0, pmeData);
+    conp.getDerivatives(numParticles, numElectrodeParticles, posData, charges, exclusions, sysToElec, elecToSys, electrodeParamArray, grad0, pme);
 
     // Project the initial gradient with preconditioning.
     if (precond) {
@@ -362,7 +356,7 @@ void ReferenceConstantPotentialCGSolver::solve(
         for (int ii = 0; ii < numElectrodeParticles; ii++) {
             charges[elecToSys[ii]] = qStep[ii];
         }
-        conp.getDerivatives(numParticles, numElectrodeParticles, posData, charges, exclusions, sysToElec, elecToSys, electrodeParamArray, gradStep, pmeData);
+        conp.getDerivatives(numParticles, numElectrodeParticles, posData, charges, exclusions, sysToElec, elecToSys, electrodeParamArray, gradStep, pme);
         for (int ii = 0; ii < numElectrodeParticles; ii++) {
             gradStep[ii] -= grad0[ii];
         }
@@ -417,7 +411,7 @@ void ReferenceConstantPotentialCGSolver::solve(
             for (int ii = 0; ii < numElectrodeParticles; ii++) {
                 charges[elecToSys[ii]] = q[ii];
             }
-            conp.getDerivatives(numParticles, numElectrodeParticles, posData, charges, exclusions, sysToElec, elecToSys, electrodeParamArray, grad, pmeData);
+            conp.getDerivatives(numParticles, numElectrodeParticles, posData, charges, exclusions, sysToElec, elecToSys, electrodeParamArray, grad, pme);
         }
         else {
             for (int ii = 0; ii < numElectrodeParticles; ii++) {
@@ -542,13 +536,9 @@ void ReferenceConstantPotential::execute(
     double* energy,
     ReferenceConstantPotentialSolver* solver
 ) {
-    pme_t pmeData;
-    pme_init(&pmeData, ewaldAlpha, numParticles, gridSize, 5, 1);
-
-    solver->solve(numParticles, numElectrodeParticles, posData, charges, exclusions, sysToElec, elecToSys, electrodeParamArray, *this, pmeData);
-    getEnergyForces(numParticles, numElectrodeParticles, posData, forceData, charges, exclusions, sysToElec, elecToSys, electrodeParamArray, energy, pmeData);
-
-    pme_destroy(pmeData);
+    ReferencePME pme(ewaldAlpha, numParticles, gridSize, 5, 1);
+    solver->solve(numParticles, numElectrodeParticles, posData, charges, exclusions, sysToElec, elecToSys, electrodeParamArray, *this, pme);
+    getEnergyForces(numParticles, numElectrodeParticles, posData, forceData, charges, exclusions, sysToElec, elecToSys, electrodeParamArray, energy, pme);
 }
 
 void ReferenceConstantPotential::getCharges(
@@ -562,12 +552,8 @@ void ReferenceConstantPotential::getCharges(
     const std::vector<std::array<double, 3> >& electrodeParamArray,
     ReferenceConstantPotentialSolver* solver
 ) {
-    pme_t pmeData;
-    pme_init(&pmeData, ewaldAlpha, numParticles, gridSize, 5, 1);
-
-    solver->solve(numParticles, numElectrodeParticles, posData, charges, exclusions, sysToElec, elecToSys, electrodeParamArray, *this, pmeData);
-
-    pme_destroy(pmeData);
+    ReferencePME pme(ewaldAlpha, numParticles, gridSize, 5, 1);
+    solver->solve(numParticles, numElectrodeParticles, posData, charges, exclusions, sysToElec, elecToSys, electrodeParamArray, *this, pme);
 }
 
 void ReferenceConstantPotential::getEnergyForces(
@@ -581,7 +567,7 @@ void ReferenceConstantPotential::getEnergyForces(
     const std::vector<int>& elecToSys,
     const std::vector<std::array<double, 3> >& electrodeParamArray,
     double* energy,
-    pme_t pmeData
+    ReferencePME& pme
 ) {
     const double SQRT_PI = sqrt(PI_M);
     const double TWO_OVER_SQRT_PI = 2.0 / SQRT_PI;
@@ -671,7 +657,7 @@ void ReferenceConstantPotential::getEnergyForces(
     
     // Reciprocal space.
     double pmeEnergy = 0.0;
-    pme_exec(pmeData, posData, forceData, charges, boxVectors, &pmeEnergy);
+    pme.exec(posData, forceData, charges, boxVectors, pmeEnergy);
     energyAccum += pmeEnergy;
 
     // Ewald self-interaction and external field contributions (all particles).
@@ -710,7 +696,7 @@ void ReferenceConstantPotential::getDerivatives(
     const std::vector<int>& elecToSys,
     const std::vector<std::array<double, 3> >& electrodeParamArray,
     std::vector<double>& chargeDerivatives,
-    pme_t pmeData
+    ReferencePME& pme
 ) {
     const double SQRT_PI = sqrt(PI_M);
     const double SELF_ALPHA_SCALE = 2.0 * ONE_4PI_EPS0 / SQRT_PI;
@@ -757,7 +743,7 @@ void ReferenceConstantPotential::getDerivatives(
     }
 
     // Reciprocal space.
-    pme_exec_charge_derivatives(pmeData, posData, chargeDerivatives, elecToSys, charges, boxVectors);
+    pme.exec_charge_derivatives(posData, chargeDerivatives, elecToSys, charges, boxVectors);
 
     // Ewald neutralizing plasma precalculation.
     double qTotal = 0.0;

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp

@@ -243,33 +243,28 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<Vec3>& a
     // PME
 
     if (pme && includeReciprocal) {
-        pme_t          pmedata; /* abstract handle for PME data */
-
-        pme_init(&pmedata,alphaEwald,numberOfAtoms,meshDim,5,1);
+        ReferencePME pme(alphaEwald,numberOfAtoms,meshDim,5,1);
 
         vector<double> charges(numberOfAtoms);
         for (int i = 0; i < numberOfAtoms; i++)
             charges[i] = atomParameters[i][QIndex];
-        pme_exec(pmedata,atomCoordinates,forces,charges,periodicBoxVectors,&recipEnergy);
+        pme.exec(atomCoordinates, forces, charges, periodicBoxVectors, recipEnergy);
 
         if (totalEnergy)
             *totalEnergy += recipEnergy;
 
-        pme_destroy(pmedata);
-
         if (ljpme) {
             // Dispersion reciprocal space terms
-            pme_init(&pmedata,alphaDispersionEwald,numberOfAtoms,dispersionMeshDim,5,1);
+            ReferencePME ljpme(alphaDispersionEwald,numberOfAtoms,dispersionMeshDim,5,1);
 
             std::vector<Vec3> dpmeforces(numberOfAtoms);
             for (int i = 0; i < numberOfAtoms; i++)
                 charges[i] = 8.0*pow(atomParameters[i][SigIndex], 3.0) * atomParameters[i][EpsIndex];
-            pme_exec_dpme(pmedata,atomCoordinates,dpmeforces,charges,periodicBoxVectors,&recipDispersionEnergy);
+            ljpme.exec_dpme(atomCoordinates, dpmeforces, charges, periodicBoxVectors, recipDispersionEnergy);
             for (int i = 0; i < numberOfAtoms; i++)
                 forces[i] += dpmeforces[i];
             if (totalEnergy)
                 *totalEnergy += recipDispersionEnergy;
-            pme_destroy(pmedata);
         }
     }
     // Ewald method

platforms/reference/src/SimTKReference/ReferencePME.cpp

@@ -29,10 +29,8 @@
  * POSSIBILITY OF SUCH DAMAGE.
  */
 
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <complex>
+#include <cmath>
+#include <vector>
 
 #include "ReferencePME.h"
 #include "ReferenceForce.h"
@@ -45,171 +43,72 @@
 
 using namespace std;
 
-typedef int    ivec[3];
-
 namespace OpenMM {
 
-struct pme
-{
-    int          natoms;
-    double       ewaldcoeff;
-
-    complex<double>* grid;             /* Memory for the grid we spread charges on.
-                                        * Element (i,j,k) is accessed as:
-                                        * grid[i*ngrid[1]*ngrid[2] + j*ngrid[2] + k]
-                                        */
-    int          ngrid[3];             /* Total grid dimensions (all data is complex!) */
-
-    int          order;                /* PME interpolation order. Almost always 4 */
-
-    /* Data for bspline interpolation, see the Essman PME paper */
-    double *     bsplines_moduli[3];   /* 3 pointers, to x/y/z bspline moduli, each of length ngrid[x/y/z]   */
-    double *     bsplines_theta[3];    /* each of x/y/z has length order*natoms */
-    double *     bsplines_dtheta[3];   /* each of x/y/z has length order*natoms */
-
-    ivec *       particleindex;        /* Array of length natoms. Each element is
-                                        * an ivec (3 ints) that specify the grid
-                                        * indices for that particular atom. Updated every step!
-                                        */
-    rvec *       particlefraction;     /* Array of length natoms. Fractional offset in the grid for
-                                        * each atom in all three dimensions.
-                                        */
-
-    /* Further explanation of index/fraction:
-     *
-     * Assume we have a cell of size 10*10*10nm, and a total grid dimension of 100*100*100 cells.
-     * In other words, each cell is 0.1*0.1*0.1 nm.
-     *
-     * If particle i has coordinates { 0.543 , 6.235 , -0.73 }, we will get:
-     *
-     * particleindex[i]    = { 5 , 62 , 92 }         (-0.73 + 10 = 9.27, we always apply PBC for grid calculations!)
-     * particlefraction[i] = { 0.43 , 0.35 , 0.7 }   (this is the fraction of the cell length where the atom is)
-     *
-     * (The reason for precaculating / storing these is that it gets a bit more complex for triclinic cells :-)
-     *
-     * In the current code version we might assume that a particle is not more than a whole box length away from
-     * the central cell, i.e., in this case we would assume all coordinates fall in -10 nm < x,y,z < 20 nm.
-     */
-
-    double       epsilon_r;             /* Dielectric coefficient to use, typically 1.0 */
-};
-
-
-/* Internal setup routines */
-
-
-
-/* Only called once from init_pme(), performance does not matter! */
-static void
-pme_calculate_bsplines_moduli(pme_t pme)
-{
-    int       nmax;
-    int       i,j,k,l,d;
-    int       order;
-    int       ndata;
-    double *  data;
-    double *  ddata;
-    double *  bsplines_data;
-    double    div;
-    double    sc,ss,arg;
-
-    nmax = 0;
-    for (d=0;d<3;d++)
-    {
-        nmax = (pme->ngrid[d] > nmax) ? pme->ngrid[d] : nmax;
-        pme->bsplines_moduli[d] = (double *) malloc(sizeof(double)*pme->ngrid[d]);
+/* Only called once from constructor, performance does not matter! */
+void ReferencePME::calculate_bsplines_moduli() {
+    int nmax = 0;
+    for (int d = 0; d < 3; d++) {
+        nmax = (ngrid[d] > nmax) ? ngrid[d] : nmax;
+        bsplines_moduli[d].resize(ngrid[d]);
     }
 
-    order = pme->order;
-
     /* temp storage in this routine */
-    data          = (double *) malloc(sizeof(double)*order);
-    ddata         = (double *) malloc(sizeof(double)*order);
-    bsplines_data = (double *) malloc(sizeof(double)*nmax);
+    vector<double> data(order);
+    vector<double> ddata(order);
+    vector<double> bsplines_data(nmax);
 
     data[order-1]=0;
     data[1]=0;
     data[0]=1;
 
-    for (k=3;k<order;k++)
-    {
-        div=1.0/(k-1.0);
+    for (int k = 3; k < order; k++) {
+        double div=1.0/(k-1.0);
         data[k-1]=0;
-        for (l=1;l<(k-1);l++)
-        {
+        for (int l = 1; l < (k-1); l++)
             data[k-l-1]=div*(l*data[k-l-2]+(k-l)*data[k-l-1]);
-        }
         data[0]=div*data[0];
     }
 
     /* differentiate */
     ddata[0]=-data[0];
-    for (k=1;k<order;k++)
-    {
+    for (int k = 1 ; k < order; k++)
         ddata[k]=data[k-1]-data[k];
-    }
 
-    div=1.0/(order-1);
+    double div=1.0/(order-1);
     data[order-1]=0;
 
-    for (l=1;l<(order-1);l++)
-    {
+    for (int l = 1; l < (order-1); l++)
         data[order-l-1]=div*(l*data[order-l-2]+(order-l)*data[order-l-1]);
-    }
     data[0]=div*data[0];
 
-    for (i=0;i<nmax;i++)
-    {
+    for (int i = 0; i < nmax; i++)
         bsplines_data[i]=0;
-    }
-    for (i=1;i<=order;i++)
-    {
+    for (int i = 1; i <= order; i++)
         bsplines_data[i]=data[i-1];
-    }
 
     /* Evaluate the actual bspline moduli for X/Y/Z */
-    for (d=0;d<3;d++)
-    {
-        ndata = pme->ngrid[d];
-        for (i=0;i<ndata;i++)
-        {
-            sc=ss=0;
-            for (j=0;j<ndata;j++)
-            {
-                arg=(2.0*M_PI*i*j)/ndata;
+    for (int d = 0; d < 3; d++) {
+        int ndata = ngrid[d];
+        for (int i = 0; i < ndata; i++) {
+            double sc = 0;
+            double ss = 0;
+            for (int j = 0; j < ndata; j++) {
+                double arg=(2.0*M_PI*i*j)/ndata;
                 sc+=bsplines_data[j]*cos(arg);
                 ss+=bsplines_data[j]*sin(arg);
             }
-            pme->bsplines_moduli[d][i]=sc*sc+ss*ss;
-        }
-        for (i=0;i<ndata;i++)
-        {
-            if (pme->bsplines_moduli[d][i]<1.0e-7)
-            {
-                pme->bsplines_moduli[d][i]=(pme->bsplines_moduli[d][(i-1+ndata)%ndata]+pme->bsplines_moduli[d][(i+1)%ndata])/2;
+            bsplines_moduli[d][i]=sc*sc+ss*ss;
         }
+        for (int i = 0; i < ndata; i++) {
+            if (bsplines_moduli[d][i]<1.0e-7)
+                bsplines_moduli[d][i]=(bsplines_moduli[d][(i-1+ndata)%ndata]+bsplines_moduli[d][(i+1)%ndata])/2;
         }
     }
-
-    /* Release temp storage */
-    free(data);
-    free(ddata);
-    free(bsplines_data);
 }
 
-static void
-pme_update_grid_index_and_fraction(pme_t    pme,
-                                   const vector<Vec3>& atomCoordinates,
-                                   const Vec3 periodicBoxVectors[3],
-                                   const Vec3 recipBoxVectors[3])
-{
-    int    i;
-    int    d;
-    double t;
-    int    ti;
-
-    for (i=0;i<pme->natoms;i++)
-    {
+void ReferencePME::update_grid_index_and_fraction(const vector<Vec3>& atomCoordinates, const Vec3 recipBoxVectors[3]) {
+    for (int i = 0; i < natoms; i++) {
         /* Index calculation (Look mom, no conditionals!):
          *
          * Both for Cuda and modern CPUs it is nice to avoid conditionals, but we still need to apply periodic boundary conditions.
@@ -247,14 +146,13 @@ pme_update_grid_index_and_fraction(pme_t    pme,
          *    (And, by adding 100.0 box lengths, we would lose a bit of numerical accuracy here!)
          */
         Vec3 coord = atomCoordinates[i];
-        for (d=0;d<3;d++)
-        {
-            t = coord[0]*recipBoxVectors[0][d]+coord[1]*recipBoxVectors[1][d]+coord[2]*recipBoxVectors[2][d];
-            t = (t-floor(t))*pme->ngrid[d];
-            ti = (int) t;
-
-            pme->particlefraction[i][d] = t - ti;
-            pme->particleindex[i][d]    = ti % pme->ngrid[d];
+        for (int d = 0; d < 3; d++) {
+            double t = coord[0]*recipBoxVectors[0][d]+coord[1]*recipBoxVectors[1][d]+coord[2]*recipBoxVectors[2][d];
+            t = (t-floor(t))*ngrid[d];
+            int ti = (int) t;
+
+            particlefraction[i][d] = t - ti;
+            particleindex[i][d]    = ti % ngrid[d];
         }
     }
 }
@@ -265,97 +163,60 @@ pme_update_grid_index_and_fraction(pme_t    pme,
  *
  * In practice, it might help to require order=4 for the cuda port.
  */
-static void
-pme_update_bsplines(pme_t    pme)
-{
-    int       i,j,k,l;
-    int       order;
-    double    dr,div;
-    double *  data;
-    double *  ddata;
-
-    order = pme->order;
-
-    for (i=0; (i<pme->natoms); i++)
-    {
-        for (j=0; j<3; j++)
-        {
+void ReferencePME::update_bsplines() {
+    for (int i = 0; i < natoms; i++) {
+        for (int j = 0; j < 3; j++) {
             /* dr is relative offset from lower cell limit */
-            dr = pme->particlefraction[i][j];
+            double dr = particlefraction[i][j];
 
-            data  = &(pme->bsplines_theta[j][i*order]);
-            ddata = &(pme->bsplines_dtheta[j][i*order]);
+            double* data  = &(bsplines_theta[j][i*order]);
+            double* ddata = &(bsplines_dtheta[j][i*order]);
             data[order-1] = 0;
             data[1]       = dr;
             data[0]       = 1-dr;
 
-            for (k=3; k<order; k++)
-            {
-                div = 1.0/(k-1.0);
+            for (int k = 3; k < order; k++) {
+                double div = 1.0/(k-1.0);
                 data[k-1] = div*dr*data[k-2];
-                for (l=1; l<(k-1); l++)
-                {
+                for (int l = 1; l < (k-1); l++)
                     data[k-l-1] = div*((dr+l)*data[k-l-2]+(k-l-dr)*data[k-l-1]);
-                }
                 data[0] = div*(1-dr)*data[0];
             }
 
             /* differentiate */
             ddata[0] = -data[0];
 
-            for (k=1; k<order; k++)
-            {
+            for (int k = 1; k < order; k++)
                 ddata[k] = data[k-1]-data[k];
-            }
 
-            div           = 1.0/(order-1);
+            double div = 1.0/(order-1);
             data[order-1] = div*dr*data[order-2];
 
-            for (l=1; l<(order-1); l++)
-            {
+            for (int l = 1; l < (order-1); l++)
                 data[order-l-1] = div*((dr+l)*data[order-l-2]+(order-l-dr)*data[order-l-1]);
-            }
             data[0] = div*(1-dr)*data[0];
         }
     }
 }
 
 
-static void
-pme_grid_spread_charge(pme_t pme, const vector<double>& charges)
-{
-    int       order;
-    int       i;
-    int       ix,iy,iz;
-    int       x0index,y0index,z0index;
-    int       xindex,yindex,zindex;
-    int       index;
-    double    q;
-    double *  thetax;
-    double *  thetay;
-    double *  thetaz;
-
-    order = pme->order;
-
+void ReferencePME::grid_spread_charge(const vector<double>& charges) {
     /* Reset the grid */
-    for (i=0;i<pme->ngrid[0]*pme->ngrid[1]*pme->ngrid[2];i++)
-    {
-        pme->grid[i] = complex<double>(0, 0);
-    }
+    for (int i = 0; i < ngrid[0]*ngrid[1]*ngrid[2]; i++)
+        grid[i] = complex<double>(0, 0);
 
-    for (i=0;i<pme->natoms;i++)
-    {
-        q = charges[i];
+    for (int i = 0; i < natoms; i++) {
+        double q = charges[i];
 
         /* Grid index for the actual atom position */
-        x0index = pme->particleindex[i][0];
-        y0index = pme->particleindex[i][1];
-        z0index = pme->particleindex[i][2];
+        int x0index = particleindex[i][0];
+        int y0index = particleindex[i][1];
+        int z0index = particleindex[i][2];
 
         /* Bspline factors for this atom in each dimension , calculated from fractional coordinates */
-        thetax  = &(pme->bsplines_theta[0][i*order]);
-        thetay  = &(pme->bsplines_theta[1][i*order]);
-        thetaz  = &(pme->bsplines_theta[2][i*order]);
+        double* thetax  = &(bsplines_theta[0][i*order]);
+        double* thetay  = &(bsplines_theta[1][i*order]);
+        double* thetaz  = &(bsplines_theta[2][i*order]);
 
         /* Loop over norder*norder*norder (typically 5*5*5) neighbor cells.
          *
@@ -375,23 +236,20 @@ pme_grid_spread_charge(pme_t pme, const vector<double>& charges)
          * 3) When we parallelize things, we only need to communicate in one direction instead of two!
          */
 
-        for (ix=0;ix<order;ix++)
-        {
+        for (int ix = 0; ix < order; ix++) {
             /* Calculate index, apply PBC so we spread to index 0/1/2 when a particle is close to the upper limit of the grid */
-            xindex = (x0index + ix) % pme->ngrid[0];
+            int xindex = (x0index + ix) % ngrid[0];
 
-            for (iy=0;iy<order;iy++)
-            {
-                yindex = (y0index + iy) % pme->ngrid[1];
+            for (int iy = 0; iy < order; iy++) {
+                int yindex = (y0index + iy) % ngrid[1];
 
-                for (iz=0;iz<order;iz++)
-                {
+                for (int iz = 0; iz < order; iz++) {
                     /* Can be optimized, but we keep it simple here */
-                    zindex            = (z0index + iz) % pme->ngrid[2];
+                    int zindex = (z0index + iz) % ngrid[2];
                     /* Calculate index in the charge grid */
-                    index             = xindex*pme->ngrid[1]*pme->ngrid[2] + yindex*pme->ngrid[2] + zindex;
+                    int index = xindex*ngrid[1]*ngrid[2] + yindex*ngrid[2] + zindex;
                     /* Add the charge times the bspline spread/interpolation factors to this grid position */
-                    pme->grid[index] += q*thetax[ix]*thetay[iy]*thetaz[iz];
+                    grid[index] += q*thetax[ix]*thetay[iy]*thetaz[iz];
                 }
             }
         }
@@ -400,67 +258,36 @@ pme_grid_spread_charge(pme_t pme, const vector<double>& charges)
 
 
 
-static void
-pme_reciprocal_convolution(pme_t     pme,
-                           const Vec3 periodicBoxVectors[3],
-                           const Vec3 recipBoxVectors[3],
-                           double *  energy)
-{
-    int kx,ky,kz;
-    int nx,ny,nz;
-    double mx,my,mz;
-    double mhx,mhy,mhz,m2;
-    double one_4pi_eps;
-    double virxx,virxy,virxz,viryy,viryz,virzz;
-    double bx,by,bz;
-    double d1,d2;
-    double eterm,vfactor,struct2,ets2;
-    double esum;
-    double factor;
-    double denom;
-    double boxfactor;
-    double maxkx,maxky,maxkz;
-
-    complex<double> *ptr;
-
-    nx = pme->ngrid[0];
-    ny = pme->ngrid[1];
-    nz = pme->ngrid[2];
-
-    one_4pi_eps = ONE_4PI_EPS0/pme->epsilon_r;
-    factor = M_PI*M_PI/(pme->ewaldcoeff*pme->ewaldcoeff);
-    boxfactor = M_PI*periodicBoxVectors[0][0]*periodicBoxVectors[1][1]*periodicBoxVectors[2][2];
-
-    esum = 0;
-    virxx = 0;
-    virxy = 0;
-    virxz = 0;
-    viryy = 0;
-    viryz = 0;
-    virzz = 0;
-
-    maxkx = (nx+1)/2;
-    maxky = (ny+1)/2;
-    maxkz = (nz+1)/2;
-
-    for (kx=0;kx<nx;kx++)
-    {
+void ReferencePME::pme_reciprocal_convolution(const Vec3 periodicBoxVectors[3], const Vec3 recipBoxVectors[3], double& energy) {
+    int nx = ngrid[0];
+    int ny = ngrid[1];
+    int nz = ngrid[2];
+
+    double one_4pi_eps = ONE_4PI_EPS0/epsilon_r;
+    double factor = M_PI*M_PI/(ewaldcoeff*ewaldcoeff);
+    double boxfactor = M_PI*periodicBoxVectors[0][0]*periodicBoxVectors[1][1]*periodicBoxVectors[2][2];
+
+    double esum = 0;
+
+    double maxkx = (nx+1)/2;
+    double maxky = (ny+1)/2;
+    double maxkz = (nz+1)/2;
+
+    for (int kx = 0; kx < nx; kx++) {
         /* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
-        mx  = (kx<maxkx) ? kx : (kx-nx);
-        mhx = mx*recipBoxVectors[0][0];
-        bx  = boxfactor*pme->bsplines_moduli[0][kx];
+        double mx  = (kx<maxkx) ? kx : (kx-nx);
+        double mhx = mx*recipBoxVectors[0][0];
+        double bx  = boxfactor*bsplines_moduli[0][kx];
 
-        for (ky=0;ky<ny;ky++)
-        {
+        for (int ky = 0; ky < ny; ky++) {
             /* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
-            my  = (ky<maxky) ? ky : (ky-ny);
-            mhy = mx*recipBoxVectors[1][0]+my*recipBoxVectors[1][1];
-            by  = pme->bsplines_moduli[1][ky];
+            double my  = (ky<maxky) ? ky : (ky-ny);
+            double mhy = mx*recipBoxVectors[1][0]+my*recipBoxVectors[1][1];
+            double by  = bsplines_moduli[1][ky];
 
-            for (kz=0;kz<nz;kz++)
-            {
+            for (int kz = 0; kz < nz; kz++) {
                 /* Pointer to the grid cell in question */
-                ptr       = pme->grid + kx*ny*nz + ky*nz + kz;
+                complex<double>& ptr = grid[kx*ny*nz + ky*nz + kz];
 
                 /* The zero frequency term is undefined due to division by the frequency below.  Set this term to zero;
                  * in the case that the net charge of the system is non-zero, this is equivalent to applying a uniform
@@ -468,228 +295,170 @@ pme_reciprocal_convolution(pme_t     pme,
                  * this neutralizing plasma is applied elsewhere.  If this term is not zeroed, however, energies and
                  * forces will be unaffected but charge derivatives for non-neutral systems will be incorrect!
                  */
-                if (kx==0 && ky==0 && kz==0)
-                {
-                    *ptr = 0;
+                if (kx==0 && ky==0 && kz==0) {
+                    ptr = 0;
                     continue;
                 }
 
                 /* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
-                mz        = (kz<maxkz) ? kz : (kz-nz);
-                mhz       = mx*recipBoxVectors[2][0]+my*recipBoxVectors[2][1]+mz*recipBoxVectors[2][2];
+                double mz        = (kz<maxkz) ? kz : (kz-nz);
+                double mhz       = mx*recipBoxVectors[2][0]+my*recipBoxVectors[2][1]+mz*recipBoxVectors[2][2];
 
                 /* Get grid data for this frequency */
-                d1        = ptr->real();
-                d2        = ptr->imag();
+                double d1        = ptr.real();
+                double d2        = ptr.imag();
 
                 /* Calculate the convolution - see the Essman/Darden paper for the equation! */
-                m2        = mhx*mhx+mhy*mhy+mhz*mhz;
-                bz        = pme->bsplines_moduli[2][kz];
-                denom     = m2*bx*by*bz;
+                double m2        = mhx*mhx+mhy*mhy+mhz*mhz;
+                double bz        = bsplines_moduli[2][kz];
+                double denom     = m2*bx*by*bz;
 
-                eterm     = one_4pi_eps*exp(-factor*m2)/denom;
+                double eterm     = one_4pi_eps*exp(-factor*m2)/denom;
 
                 /* write back convolution data to grid */
-                ptr->real(d1*eterm);
-                ptr->imag(d2*eterm);
+                ptr.real(d1*eterm);
+                ptr.imag(d2*eterm);
 
-                struct2   = (d1*d1+d2*d2);
+                double struct2   = (d1*d1+d2*d2);
 
                 /* Long-range PME contribution to the energy for this frequency */
-                ets2      = eterm*struct2;
+                double ets2      = eterm*struct2;
                 esum     += ets2;
             }
         }
     }
 
     /* The factor 0.5 is nothing special, but it is better to have it here than inside the loop :-) */
-    *energy = 0.5*esum;
+    energy = 0.5*esum;
 }
 
 
-static void
-dpme_reciprocal_convolution(pme_t pme,
-                           const Vec3 periodicBoxVectors[3],
-                           const Vec3 recipBoxVectors[3],
-                           double* energy)
-{
-    int kx,ky,kz;
-    int nx,ny,nz;
-    double mx,my,mz;
-    double mhx,mhy,mhz,m2;
-    double bx,by,bz;
-    double d1,d2;
-    double eterm,struct2,ets2;
-    double esum;
-    double denom;
-    double boxfactor;
-    double maxkx,maxky,maxkz;
+void ReferencePME::dpme_reciprocal_convolution(const Vec3 periodicBoxVectors[3], const Vec3 recipBoxVectors[3], double& energy) {
+    int nx = ngrid[0];
+    int ny = ngrid[1];
+    int nz = ngrid[2];
 
-    complex<double> *ptr;
+    double boxfactor = -2*M_PI*sqrt(M_PI) / (6.0*periodicBoxVectors[0][0]*periodicBoxVectors[1][1]*periodicBoxVectors[2][2]);
 
-    nx = pme->ngrid[0];
-    ny = pme->ngrid[1];
-    nz = pme->ngrid[2];
+    double esum = 0;
 
-    boxfactor = -2*M_PI*sqrt(M_PI) / (6.0*periodicBoxVectors[0][0]*periodicBoxVectors[1][1]*periodicBoxVectors[2][2]);
+    double maxkx = (nx+1)/2;
+    double maxky = (ny+1)/2;
+    double maxkz = (nz+1)/2;
 
-    esum = 0;
-
-    maxkx = (nx+1)/2;
-    maxky = (ny+1)/2;
-    maxkz = (nz+1)/2;
-
-    double bfac = M_PI / pme->ewaldcoeff;
+    double bfac = M_PI / ewaldcoeff;
     double fac1 = 2.0*M_PI*M_PI*M_PI*sqrt(M_PI);
-    double fac2 = pme->ewaldcoeff*pme->ewaldcoeff*pme->ewaldcoeff;
-    double fac3 = -2.0*pme->ewaldcoeff*M_PI*M_PI;
-    double b, m, m3, expfac, expterm, erfcterm;
+    double fac2 = ewaldcoeff*ewaldcoeff*ewaldcoeff;
+    double fac3 = -2.0*ewaldcoeff*M_PI*M_PI;
 
-    for (kx=0;kx<nx;kx++)
-    {
+    for (int kx = 0; kx < nx; kx++) {
         /* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
-        mx  = ((kx<maxkx) ? kx : (kx-nx));
-        mhx = mx*recipBoxVectors[0][0];
-        bx  = pme->bsplines_moduli[0][kx];
+        double mx  = ((kx<maxkx) ? kx : (kx-nx));
+        double mhx = mx*recipBoxVectors[0][0];
+        double bx  = bsplines_moduli[0][kx];
 
-        for (ky=0;ky<ny;ky++)
-        {
+        for (int ky = 0; ky < ny; ky++) {
             /* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
-            my  = ((ky<maxky) ? ky : (ky-ny));
-            mhy = mx*recipBoxVectors[1][0]+my*recipBoxVectors[1][1];
-            by  = pme->bsplines_moduli[1][ky];
+            double my  = ((ky<maxky) ? ky : (ky-ny));
+            double mhy = mx*recipBoxVectors[1][0]+my*recipBoxVectors[1][1];
+            double by  = bsplines_moduli[1][ky];
 
-            for (kz=0;kz<nz;kz++)
-            {
+            for (int kz = 0; kz < nz; kz++) {
                 /*
                  * Unlike the Coulombic case, there's an m=0 term so all terms are considered here.
                  */
 
                 /* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
-                mz        = ((kz<maxkz) ? kz : (kz-nz));
-                mhz       = mx*recipBoxVectors[2][0]+my*recipBoxVectors[2][1]+mz*recipBoxVectors[2][2];
+                double mz  = ((kz<maxkz) ? kz : (kz-nz));
+                double mhz = mx*recipBoxVectors[2][0]+my*recipBoxVectors[2][1]+mz*recipBoxVectors[2][2];
 
                 /* Pointer to the grid cell in question */
-                ptr       = pme->grid + kx*ny*nz + ky*nz + kz;
+                complex<double>& ptr = grid[kx*ny*nz + ky*nz + kz];
 
                 /* Get grid data for this frequency */
-                d1        = ptr->real();
-                d2        = ptr->imag();
+                double d1 = ptr.real();
+                double d2 = ptr.imag();
 
                 /* Calculate the convolution - see the Essman/Darden paper for the equation! */
-                m2        = mhx*mhx+mhy*mhy+mhz*mhz;
-                bz        = pme->bsplines_moduli[2][kz];
-                denom     = boxfactor / (bx*by*bz);
+                double m2    = mhx*mhx+mhy*mhy+mhz*mhz;
+                double bz    = bsplines_moduli[2][kz];
+                double denom = boxfactor / (bx*by*bz);
 
-                m = sqrt(m2);
-                m3 = m*m2;
-                b = bfac*m;
-                expfac = -b*b;
-                erfcterm = erfc(b);
-                expterm = exp(expfac);
+                double m = sqrt(m2);
+                double m3 = m*m2;
+                double b = bfac*m;
+                double expfac = -b*b;
+                double erfcterm = erfc(b);
+                double expterm = exp(expfac);
 
-                eterm     = (fac1*erfcterm*m3 + expterm*(fac2 + fac3*m2)) * denom;
+                double eterm = (fac1*erfcterm*m3 + expterm*(fac2 + fac3*m2)) * denom;
 
                 /* write back convolution data to grid */
-                ptr->real(d1*eterm);
-                ptr->imag(d2*eterm);
+                ptr.real(d1*eterm);
+                ptr.imag(d2*eterm);
 
-                struct2   = (d1*d1+d2*d2);
+                double struct2 = (d1*d1+d2*d2);
 
                 /* Long-range PME contribution to the energy for this frequency */
-                ets2      = eterm*struct2;
+                double ets2 = eterm*struct2;
                 esum += ets2;
             }
         }
     }
     // Remember the C6 energy is attractive, hence the negative sign.
-    *energy = 0.5*esum;
+    energy = 0.5*esum;
 }
 
 
-static void
-pme_grid_interpolate_force(pme_t pme,
-                           const Vec3 recipBoxVectors[3],
-                           const vector<double>& charges,
-                           vector<Vec3>& forces)
-{
-    int       i;
-    int       ix,iy,iz;
-    int       x0index,y0index,z0index;
-    int       xindex,yindex,zindex;
-    int       index;
-    int       order;
-    double    q;
-    double *  thetax;
-    double *  thetay;
-    double *  thetaz;
-    double *  dthetax;
-    double *  dthetay;
-    double *  dthetaz;
-    double    tx,ty,tz;
-    double    dtx,dty,dtz;
-    double    fx,fy,fz;
-    double    gridvalue;
-    int       nx,ny,nz;
-
-    nx    = pme->ngrid[0];
-    ny    = pme->ngrid[1];
-    nz    = pme->ngrid[2];
-
-    order = pme->order;
-
+void ReferencePME::grid_interpolate_force(const Vec3 recipBoxVectors[3], const vector<double>& charges, vector<Vec3>& forces) {
     /* This is almost identical to the charge spreading routine! */
 
-    for (i=0;i<pme->natoms;i++)
-    {
-        fx = fy = fz = 0;
+    for (int i = 0; i < natoms; i++) {
+        double fx = 0, fy = 0, fz = 0;
 
-        q = charges[i];
+        double q = charges[i];
 
         /* Grid index for the actual atom position */
-        x0index = pme->particleindex[i][0];
-        y0index = pme->particleindex[i][1];
-        z0index = pme->particleindex[i][2];
+        int x0index = particleindex[i][0];
+        int y0index = particleindex[i][1];
+        int z0index = particleindex[i][2];
 
         /* Bspline factors for this atom in each dimension , calculated from fractional coordinates */
-        thetax  = &(pme->bsplines_theta[0][i*order]);
-        thetay  = &(pme->bsplines_theta[1][i*order]);
-        thetaz  = &(pme->bsplines_theta[2][i*order]);
-        dthetax = &(pme->bsplines_dtheta[0][i*order]);
-        dthetay = &(pme->bsplines_dtheta[1][i*order]);
-        dthetaz = &(pme->bsplines_dtheta[2][i*order]);
+        double* thetax  = &(bsplines_theta[0][i*order]);
+        double* thetay  = &(bsplines_theta[1][i*order]);
+        double* thetaz  = &(bsplines_theta[2][i*order]);
+        double* dthetax = &(bsplines_dtheta[0][i*order]);
+        double* dthetay = &(bsplines_dtheta[1][i*order]);
+        double* dthetaz = &(bsplines_dtheta[2][i*order]);
 
-        /* See pme_grid_spread_charge() for comments about the order here, and only interpolation in one direction */
+        /* See grid_spread_charge() for comments about the order here, and only interpolation in one direction */
 
         /* Since we will add order^3 (typically 5*5*5=125) terms to the force on each particle, we use temporary fx/fy/fz
          * variables, and only add it to memory forces[] at the end.
          */
-        for (ix=0;ix<order;ix++)
-        {
-            xindex = (x0index + ix) % pme->ngrid[0];
+        for (int ix = 0; ix < order; ix++) {
+            int xindex = (x0index + ix) % ngrid[0];
             /* Get both the bspline factor and its derivative with respect to the x coordinate! */
-            tx     = thetax[ix];
-            dtx    = dthetax[ix];
+            double tx  = thetax[ix];
+            double dtx = dthetax[ix];
 
-            for (iy=0;iy<order;iy++)
-            {
-                yindex = (y0index + iy) % pme->ngrid[1];
+            for (int iy = 0; iy < order; iy++) {
+                int yindex = (y0index + iy) % ngrid[1];
                 /* bspline + derivative wrt y */
-                ty     = thetay[iy];
-                dty    = dthetay[iy];
+                double ty  = thetay[iy];
+                double dty = dthetay[iy];
 
-                for (iz=0;iz<order;iz++)
-                {
+                for (int iz = 0; iz < order; iz++) {
                     /* Can be optimized, but we keep it simple here */
-                    zindex               = (z0index + iz) % pme->ngrid[2];
+                    int zindex = (z0index + iz) % ngrid[2];
                     /* bspline + derivative wrt z */
-                    tz                   = thetaz[iz];
-                    dtz                  = dthetaz[iz];
-                    index                = xindex*pme->ngrid[1]*pme->ngrid[2] + yindex*pme->ngrid[2] + zindex;
+                    double tz  = thetaz[iz];
+                    double dtz = dthetaz[iz];
+                    int index  = xindex*ngrid[1]*ngrid[2] + yindex*ngrid[2] + zindex;
 
                     /* Get the fft+convoluted+ifft:d data from the grid, which must be real by definition */
                     /* Checking that the imaginary part is indeed zero might be a good check :-) */
-                    gridvalue            = pme->grid[index].real();
+                    double gridvalue = grid[index].real();
 
                     /* The d component of the force is calculated by taking the derived bspline in dimension d, normal bsplines in the other two */
                     fx += dtx*ty*tz*gridvalue;
@@ -699,90 +468,58 @@ pme_grid_interpolate_force(pme_t pme,
             }
         }
         /* Update memory force, note that we multiply by charge and some box stuff */
-        forces[i][0] -= q*(fx*nx*recipBoxVectors[0][0]);
-        forces[i][1] -= q*(fx*nx*recipBoxVectors[1][0]+fy*ny*recipBoxVectors[1][1]);
-        forces[i][2] -= q*(fx*nx*recipBoxVectors[2][0]+fy*ny*recipBoxVectors[2][1]+fz*nz*recipBoxVectors[2][2]);
+        forces[i][0] -= q*(fx*ngrid[0]*recipBoxVectors[0][0]);
+        forces[i][1] -= q*(fx*ngrid[0]*recipBoxVectors[1][0]+fy*ngrid[1]*recipBoxVectors[1][1]);
+        forces[i][2] -= q*(fx*ngrid[0]*recipBoxVectors[2][0]+fy*ngrid[1]*recipBoxVectors[2][1]+fz*ngrid[2]*recipBoxVectors[2][2]);
     }
 }
 
 
-static void
-pme_grid_interpolate_charge_derivatives(pme_t pme,
-                                        const Vec3 recipBoxVectors[3],
-                                        const vector<double>& charges,
-                                        vector<double>& chargeDerivatives,
-                                        const vector<int>& chargeIndices)
-{
-    int       nderiv, ideriv, i;
-    int       ix,iy,iz;
-    int       x0index,y0index,z0index;
-    int       xindex,yindex,zindex;
-    int       index;
-    int       order;
-    double    q;
-    double *  thetax;
-    double *  thetay;
-    double *  thetaz;
-    double    tx,ty,tz;
-    double    dq;
-    double    gridvalue;
-    int       nx,ny,nz;
-
-    nx    = pme->ngrid[0];
-    ny    = pme->ngrid[1];
-    nz    = pme->ngrid[2];
-
-    order = pme->order;
-
-    /* This is similar to pme_grid_interpolate_force() */
-    
-    nderiv = chargeIndices.size();
-    for (ideriv=0;ideriv<nderiv;ideriv++)
-    {
-        i = chargeIndices[ideriv];
-        dq = 0;
-
-        q = charges[i];
+void ReferencePME::grid_interpolate_charge_derivatives(const Vec3 recipBoxVectors[3], const vector<double>& charges,
+            vector<double>& chargeDerivatives, const vector<int>& chargeIndices) {
+    /* This is similar to grid_interpolate_force() */
+    
+    int nderiv = chargeIndices.size();
+    for (int ideriv = 0; ideriv < nderiv; ideriv++) {
+        int i = chargeIndices[ideriv];
+        double dq = 0;
 
         /* Grid index for the actual atom position */
-        x0index = pme->particleindex[i][0];
-        y0index = pme->particleindex[i][1];
-        z0index = pme->particleindex[i][2];
+        int x0index = particleindex[i][0];
+        int y0index = particleindex[i][1];
+        int z0index = particleindex[i][2];
 
         /* Bspline factors for this atom in each dimension , calculated from fractional coordinates */
-        thetax  = &(pme->bsplines_theta[0][i*order]);
-        thetay  = &(pme->bsplines_theta[1][i*order]);
-        thetaz  = &(pme->bsplines_theta[2][i*order]);
+        double* thetax  = &(bsplines_theta[0][i*order]);
+        double* thetay  = &(bsplines_theta[1][i*order]);
+        double* thetaz  = &(bsplines_theta[2][i*order]);
 
-        /* See pme_grid_spread_charge() for comments about the order here, and only interpolation in one direction */
+        /* See grid_spread_charge() for comments about the order here, and only interpolation in one direction */
 
         /* Since we will add order^3 (typically 5*5*5=125) terms to the charge
          * derivative on each particle, we use a temporary dq variable, and only
          * add it to memory forces[] at the end.
          */
-        for (ix=0;ix<order;ix++)
-        {
-            xindex = (x0index + ix) % pme->ngrid[0];
+        for (int ix = 0; ix < order; ix++) {
+            int xindex = (x0index + ix) % ngrid[0];
             /* Get the bspline factor with respect to the x coordinate */
-            tx     = thetax[ix];
+            double tx  = thetax[ix];
 
-            for (iy=0;iy<order;iy++)
-            {
-                yindex = (y0index + iy) % pme->ngrid[1];
+            for (int iy = 0; iy < order; iy++) {
+                int yindex = (y0index + iy) % ngrid[1];
                 /* bspline wrt y */
-                ty     = thetay[iy];
+                double ty  = thetay[iy];
 
-                for (iz=0;iz<order;iz++)
-                {
+                for (int iz = 0; iz < order; iz++) {
                     /* Can be optimized, but we keep it simple here */
-                    zindex    = (z0index + iz) % pme->ngrid[2];
+                    int zindex = (z0index + iz) % ngrid[2];
                     /* bspline wrt z */
-                    tz        = thetaz[iz];
-                    index     = xindex*pme->ngrid[1]*pme->ngrid[2] + yindex*pme->ngrid[2] + zindex;
+                    double tz  = thetaz[iz];
+                    int index  = xindex*ngrid[1]*ngrid[2] + yindex*ngrid[2] + zindex;
 
                     /* Get the fft+convoluted+ifft:d data from the grid, which must be real by definition */
                     /* Checking that the imaginary part is indeed zero might be a good check :-) */
-                    gridvalue = pme->grid[index].real();
+                    double gridvalue = grid[index].real();
 
                     /* The d component of the force is calculated by taking the derived bspline in dimension d, normal bsplines in the other two */
                     dq += tx*ty*tz*gridvalue;
@@ -794,222 +531,146 @@ pme_grid_interpolate_charge_derivatives(pme_t pme,
     }
 }
 
-
-
-/* EXPORTED ROUTINES */
-
-int
-pme_init(pme_t *       ppme,
-         double        ewaldcoeff,
-         int           natoms,
-         const int     ngrid[3],
-         int           pme_order,
-         double        epsilon_r)
-{
-    pme_t pme;
-    int   d;
-
-    pme = (pme_t) malloc(sizeof(struct pme));
-
-    pme->order       = pme_order;
-    pme->epsilon_r   = epsilon_r;
-    pme->ewaldcoeff  = ewaldcoeff;
-    pme->natoms      = natoms;
-
-    for (d=0;d<3;d++)
-    {
-        pme->ngrid[d]            = ngrid[d];
-        pme->bsplines_theta[d]   = (double *)malloc(sizeof(double)*pme_order*natoms);
-        pme->bsplines_dtheta[d]  = (double *)malloc(sizeof(double)*pme_order*natoms);
+ReferencePME::ReferencePME(double ewaldcoeff, int natoms, const int ngrid[3], int pme_order, double epsilon_r) :
+                    order(pme_order), epsilon_r(epsilon_r), ewaldcoeff(ewaldcoeff), natoms(natoms) {
+    for (int d = 0; d < 3; d++) {
+        this->ngrid[d] = ngrid[d];
+        bsplines_theta[d].resize(pme_order*natoms);
+        bsplines_dtheta[d].resize(pme_order*natoms);
     }
 
-    pme->particlefraction = (rvec *)malloc(sizeof(rvec)*natoms);
-    pme->particleindex    = (ivec *)malloc(sizeof(ivec)*natoms);
+    particlefraction.resize(natoms);
+    particleindex.resize(natoms);
 
     /* Allocate charge grid storage */
-    pme->grid        = (complex<double> *)malloc(sizeof(complex<double>)*ngrid[0]*ngrid[1]*ngrid[2]);
+    grid.resize(ngrid[0]*ngrid[1]*ngrid[2]);
 
     /* Setup bspline moduli (see Essman paper) */
-    pme_calculate_bsplines_moduli(pme);
-
-    *ppme = pme;
-
-    return 0;
+    calculate_bsplines_moduli();
 }
 
-
-
-
-
-int pme_exec(pme_t       pme,
-             const vector<Vec3>& atomCoordinates,
-             vector<Vec3>& forces,
-             const vector<double>& charges,
-             const Vec3 periodicBoxVectors[3],
-             double* energy)
-{
+void ReferencePME::exec(const vector<Vec3>& atomCoordinates, vector<Vec3>& forces, const vector<double>& charges,
+            const Vec3 periodicBoxVectors[3], double& energy) {
     /* Routine is called with coordinates in x, a box, and charges in q */
 
     Vec3 recipBoxVectors[3];
     ReferenceForce::invertBoxVectors(periodicBoxVectors, recipBoxVectors);
     
     /* Before we can do the actual interpolation, we need to recalculate and update
-     * the indices for each particle in the charge grid (initialized in pme_init()),
+     * the indices for each particle in the charge grid (initialized in constructor),
      * and what its fractional offset in this grid cell is.
      */
 
     /* Update charge grid indices and fractional offsets for each atom.
      * The indices/fractions are stored internally in the pme datatype
      */
-    pme_update_grid_index_and_fraction(pme,atomCoordinates,periodicBoxVectors,recipBoxVectors);
+    update_grid_index_and_fraction(atomCoordinates, recipBoxVectors);
 
     /* Calculate bsplines (and their differentials) from current fractional coordinates, store in pme structure */
-    pme_update_bsplines(pme);
+    update_bsplines();
 
     /* Spread the charges on grid (using newly calculated bsplines in the pme structure) */
-    pme_grid_spread_charge(pme, charges);
+    grid_spread_charge(charges);
 
     /* do 3d-fft */
-    vector<size_t> shape = {(size_t) pme->ngrid[0], (size_t) pme->ngrid[1], (size_t) pme->ngrid[2]};
+    vector<size_t> shape = {(size_t) ngrid[0], (size_t) ngrid[1], (size_t) ngrid[2]};
     vector<size_t> axes = {0, 1, 2};
-    vector<ptrdiff_t> stride = {(ptrdiff_t) (pme->ngrid[1]*pme->ngrid[2]*sizeof(complex<double>)),
-                                (ptrdiff_t) (pme->ngrid[2]*sizeof(complex<double>)),
+    vector<ptrdiff_t> stride = {(ptrdiff_t) (ngrid[1]*ngrid[2]*sizeof(complex<double>)),
+                                (ptrdiff_t) (ngrid[2]*sizeof(complex<double>)),
                                 (ptrdiff_t) sizeof(complex<double>)};
-    pocketfft::c2c(shape, stride, stride, axes, true, pme->grid, pme->grid, 1.0, 0);
+    pocketfft::c2c(shape, stride, stride, axes, true, grid.data(), grid.data(), 1.0, 0);
 
     /* solve in k-space */
-    pme_reciprocal_convolution(pme,periodicBoxVectors,recipBoxVectors,energy);
+    pme_reciprocal_convolution(periodicBoxVectors, recipBoxVectors, energy);
 
     /* do 3d-invfft */
-    pocketfft::c2c(shape, stride, stride, axes, false, pme->grid, pme->grid, 1.0, 0);
+    pocketfft::c2c(shape, stride, stride, axes, false, grid.data(), grid.data(), 1.0, 0);
 
     /* Get the particle forces from the grid and bsplines in the pme structure */
-    pme_grid_interpolate_force(pme,recipBoxVectors,charges,forces);
-
-    return 0;
+    grid_interpolate_force(recipBoxVectors, charges, forces);
 }
 
-
-int pme_exec_charge_derivatives(pme_t pme,
-                                const vector<Vec3>& atomCoordinates,
-                                vector<double>& chargeDerivatives,
-                                const vector<int>& chargeIndices,
-                                const vector<double>& charges,
-                                const Vec3 periodicBoxVectors[3])
-{
+void ReferencePME::exec_charge_derivatives(const vector<Vec3>& atomCoordinates, vector<double>& chargeDerivatives,
+            const vector<int>& chargeIndices, const vector<double>& charges, const Vec3 periodicBoxVectors[3]) {
     /* Routine is called with coordinates in x, a box, and charges in q */
 
     Vec3 recipBoxVectors[3];
     ReferenceForce::invertBoxVectors(periodicBoxVectors, recipBoxVectors);
     
     /* Before we can do the actual interpolation, we need to recalculate and update
-     * the indices for each particle in the charge grid (initialized in pme_init()),
+     * the indices for each particle in the charge grid (initialized in constructor),
      * and what its fractional offset in this grid cell is.
      */
 
     /* Update charge grid indices and fractional offsets for each atom.
      * The indices/fractions are stored internally in the pme datatype
      */
-    pme_update_grid_index_and_fraction(pme,atomCoordinates,periodicBoxVectors,recipBoxVectors);
+    update_grid_index_and_fraction(atomCoordinates, recipBoxVectors);
 
     /* Calculate bsplines (and their differentials) from current fractional coordinates, store in pme structure */
-    pme_update_bsplines(pme);
+    update_bsplines();
 
     /* Spread the charges on grid (using newly calculated bsplines in the pme structure) */
-    pme_grid_spread_charge(pme, charges);
+    grid_spread_charge(charges);
 
     /* do 3d-fft */
-    vector<size_t> shape = {(size_t) pme->ngrid[0], (size_t) pme->ngrid[1], (size_t) pme->ngrid[2]};
+    vector<size_t> shape = {(size_t) ngrid[0], (size_t) ngrid[1], (size_t) ngrid[2]};
     vector<size_t> axes = {0, 1, 2};
-    vector<ptrdiff_t> stride = {(ptrdiff_t) (pme->ngrid[1]*pme->ngrid[2]*sizeof(complex<double>)),
-                                (ptrdiff_t) (pme->ngrid[2]*sizeof(complex<double>)),
+    vector<ptrdiff_t> stride = {(ptrdiff_t) (ngrid[1]*ngrid[2]*sizeof(complex<double>)),
+                                (ptrdiff_t) (ngrid[2]*sizeof(complex<double>)),
                                 (ptrdiff_t) sizeof(complex<double>)};
-    pocketfft::c2c(shape, stride, stride, axes, true, pme->grid, pme->grid, 1.0, 0);
+    pocketfft::c2c(shape, stride, stride, axes, true, grid.data(), grid.data(), 1.0, 0);
 
     /* solve in k-space */
     double energy;
-    pme_reciprocal_convolution(pme,periodicBoxVectors,recipBoxVectors,&energy);
+    pme_reciprocal_convolution(periodicBoxVectors, recipBoxVectors, energy);
 
     /* do 3d-invfft */
-    pocketfft::c2c(shape, stride, stride, axes, false, pme->grid, pme->grid, 1.0, 0);
+    pocketfft::c2c(shape, stride, stride, axes, false, grid.data(), grid.data(), 1.0, 0);
 
     /* Get the charge derivatives from the grid and bsplines in the pme structure */
-    pme_grid_interpolate_charge_derivatives(pme,recipBoxVectors,charges,chargeDerivatives,chargeIndices);
-
-    return 0;
+    grid_interpolate_charge_derivatives(recipBoxVectors, charges, chargeDerivatives, chargeIndices);
 }
 
-
-int pme_exec_dpme(pme_t       pme,
-             const vector<Vec3>& atomCoordinates,
-             vector<Vec3>& forces,
-             const vector<double>& c6s,
-             const Vec3 periodicBoxVectors[3],
-             double* energy)
-{
+void ReferencePME::exec_dpme(const vector<Vec3>& atomCoordinates, vector<Vec3>& forces, const vector<double>& c6s,
+            const Vec3 periodicBoxVectors[3], double& energy) {
     /* Routine is called with coordinates in x, a box, and charges in q */
 
     Vec3 recipBoxVectors[3];
     ReferenceForce::invertBoxVectors(periodicBoxVectors, recipBoxVectors);
 
     /* Before we can do the actual interpolation, we need to recalculate and update
-     * the indices for each particle in the charge grid (initialized in pme_init()),
+     * the indices for each particle in the charge grid (initialized in constructor),
      * and what its fractional offset in this grid cell is.
      */
 
     /* Update charge grid indices and fractional offsets for each atom.
      * The indices/fractions are stored internally in the pme datatype
      */
-    pme_update_grid_index_and_fraction(pme,atomCoordinates,periodicBoxVectors,recipBoxVectors);
+    update_grid_index_and_fraction(atomCoordinates, recipBoxVectors);
 
     /* Calculate bsplines (and their differentials) from current fractional coordinates, store in pme structure */
-    pme_update_bsplines(pme);
+    update_bsplines();
 
     /* Spread the charges on grid (using newly calculated bsplines in the pme structure) */
-    pme_grid_spread_charge(pme, c6s);
+    grid_spread_charge(c6s);
 
     /* do 3d-fft */
-    vector<size_t> shape = {(size_t) pme->ngrid[0], (size_t) pme->ngrid[1], (size_t) pme->ngrid[2]};
+    vector<size_t> shape = {(size_t) ngrid[0], (size_t) ngrid[1], (size_t) ngrid[2]};
     vector<size_t> axes = {0, 1, 2};
-    vector<ptrdiff_t> stride = {(ptrdiff_t) (pme->ngrid[1]*pme->ngrid[2]*sizeof(complex<double>)),
-                                (ptrdiff_t) (pme->ngrid[2]*sizeof(complex<double>)),
+    vector<ptrdiff_t> stride = {(ptrdiff_t) (ngrid[1]*ngrid[2]*sizeof(complex<double>)),
+                                (ptrdiff_t) (ngrid[2]*sizeof(complex<double>)),
                                 (ptrdiff_t) sizeof(complex<double>)};
-    pocketfft::c2c(shape, stride, stride, axes, true, pme->grid, pme->grid, 1.0, 0);
+    pocketfft::c2c(shape, stride, stride, axes, true, grid.data(), grid.data(), 1.0, 0);
 
     /* solve in k-space */
-    dpme_reciprocal_convolution(pme,periodicBoxVectors,recipBoxVectors,energy);
+    dpme_reciprocal_convolution(periodicBoxVectors, recipBoxVectors, energy);
 
     /* do 3d-invfft */
-    pocketfft::c2c(shape, stride, stride, axes, false, pme->grid, pme->grid, 1.0, 0);
+    pocketfft::c2c(shape, stride, stride, axes, false, grid.data(), grid.data(), 1.0, 0);
 
     /* Get the particle forces from the grid and bsplines in the pme structure */
-    pme_grid_interpolate_force(pme,recipBoxVectors,c6s,forces);
-
-    return 0;
-}
-
-
-int
-pme_destroy(pme_t    pme)
-{
-    int d;
-
-    free(pme->grid);
-
-    for (d=0;d<3;d++)
-    {
-        free(pme->bsplines_moduli[d]);
-        free(pme->bsplines_theta[d]);
-        free(pme->bsplines_dtheta[d]);
-    }
-
-    free(pme->particlefraction);
-    free(pme->particleindex);
-
-    /* destroy structure itself */
-    free(pme);
-
-    return 0;
+    grid_interpolate_force(recipBoxVectors, c6s, forces);
 }
 
 } // namespace OpenMM

platforms/reference/tests/TestReferenceEwald.cpp

@@ -59,14 +59,13 @@ void testReferencePmeDerivatives() {
     double ewaldAlpha = 1.752;
     int gridSize[3] = {54, 49, 43};
 
-    pme_t pme;
-    pme_init(&pme, ewaldAlpha, numParticles, gridSize, 5, 1);
+    ReferencePME pme(ewaldAlpha, numParticles, gridSize, 5, 1);
 
     double dummyEnergy=0;
     vector<Vec3> dummyForces(numParticles);
 
     vector<Vec3> testForces(numParticles);
-    pme_exec(pme, positions, testForces, charges, boxVectors, &dummyEnergy);
+    pme.exec(positions, testForces, charges, boxVectors, dummyEnergy);
 
     for (int i = 0; i < numParticles; i++) {
         for (int j = 0; j < 3; j++) {
@@ -74,11 +73,11 @@ void testReferencePmeDerivatives() {
 
             double energyLess = 0.0;
             positions[i][j] = referencePosition - DELTA;
-            pme_exec(pme, positions, dummyForces, charges, boxVectors, &energyLess);
+            pme.exec(positions, dummyForces, charges, boxVectors, energyLess);
 
             double energyMore = 0.0;
             positions[i][j] = referencePosition + DELTA;
-            pme_exec(pme, positions, dummyForces, charges, boxVectors, &energyMore);
+            pme.exec(positions, dummyForces, charges, boxVectors, energyMore);
 
             positions[i][j] = referencePosition;
 
@@ -87,25 +86,23 @@ void testReferencePmeDerivatives() {
     }
 
     vector<double> testDerivatives(numParticles);
-    pme_exec_charge_derivatives(pme, positions, testDerivatives, indices, charges, boxVectors);
+    pme.exec_charge_derivatives(positions, testDerivatives, indices, charges, boxVectors);
 
     for (int i = 0; i < numParticles; i++) {
         double referenceCharge = charges[i];
 
         double energyLess = 0.0;
         charges[i] = referenceCharge - DELTA;
-        pme_exec(pme, positions, dummyForces, charges, boxVectors, &energyLess);
+        pme.exec(positions, dummyForces, charges, boxVectors, energyLess);
 
         double energyMore = 0.0;
         charges[i] = referenceCharge + DELTA;
-        pme_exec(pme, positions, dummyForces, charges, boxVectors, &energyMore);
+        pme.exec(positions, dummyForces, charges, boxVectors, energyMore);
 
         charges[i] = referenceCharge;
 
         ASSERT_EQUAL_TOL((energyMore - energyLess) / (2 * DELTA), testDerivatives[i], EPSILON);
     }
-
-    pme_destroy(pme);
 }
 
 void runPlatformTests() {

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHippoNonbondedForce.cpp

@@ -2873,8 +2873,7 @@ double AmoebaReferencePmeHippoNonbondedForce::calculateDispersionPairIxn(const M
 }
 
 double AmoebaReferencePmeHippoNonbondedForce::computeReciprocalSpaceDispersionForceAndEnergy(const vector<MultipoleParticleData>& particleData, vector<Vec3>& forces) const {
-    pme_t pmedata;
-    pme_init(&pmedata, _dalphaEwald, _numParticles, _dpmeGridDimensions, 5, 1);
+    ReferencePME pme(_dalphaEwald, _numParticles, _dpmeGridDimensions, 5, 1);
     vector<double> charges(_numParticles);
     vector<Vec3> dpmeforces(_numParticles, Vec3()), coords;
     for (int i = 0; i < _numParticles; i++) {
@@ -2882,8 +2881,7 @@ double AmoebaReferencePmeHippoNonbondedForce::computeReciprocalSpaceDispersionFo
         coords.push_back(particleData[i].position);
     }
     double recipDispersionEnergy;
-    pme_exec_dpme(pmedata, coords, dpmeforces, charges, _periodicBoxVectors, &recipDispersionEnergy);
-    pme_destroy(pmedata);
+    pme.exec_dpme(coords, dpmeforces, charges, _periodicBoxVectors, recipDispersionEnergy);
     for (int i = 0; i < _numParticles; i++)
         forces[i] += dpmeforces[i];
     return recipDispersionEnergy;

plugins/cpupme/tests/TestCpuPme.cpp

@@ -650,12 +650,10 @@ void testPME(bool triclinic, bool nonNeutral) {
 
     // Get charge derivatives from the reference PME implementation.
 
-    pme_t referencePme;
     int gridSize[3] = {gridx, gridy, gridz};
-    pme_init(&referencePme, alpha, numParticles, gridSize, 5, 1);
+    ReferencePME referencePme(alpha, numParticles, gridSize, 5, 1);
     vector<double> testDerivatives(numParticles);
-    pme_exec_charge_derivatives(referencePme, positions, testDerivatives, testIndices, testCharges, boxVectors);
-    pme_destroy(referencePme);
+    referencePme.exec_charge_derivatives(positions, testDerivatives, testIndices, testCharges, boxVectors);
 
     // See if they match.
     for (int i = 0; i < numParticles; i++) {