C++ interface for surface area calculation
authorTeemu Murtola <teemu.murtola@gmail.com>
Sat, 29 Nov 2014 04:58:51 +0000 (06:58 +0200)
committerTeemu Murtola <teemu.murtola@gmail.com>
Sun, 14 Dec 2014 15:37:01 +0000 (16:37 +0100)
Convert nsc_dclm_pbc() to a C++ class that provides the same interface,
except that the number of surface dots is specified separately, not for
each calculation.  The internal implementation is still mainly the same,
but the new interface makes it possible to get rid of the global
variables storing the unit sphere surface dots.  This in turn is
requires for making the algorithm thread-safe.

Clean up the implementation a bit by making internal functions static
(they will later become members of the Impl class), and replacing some
error-handling-by-writing-to-stderr-using-global-variables with asserts
or silent proper behavior.

Add tests for all the different paths of the unit sphere surface dot
generation.  These help ensure that the algorithm stays unmodified when
it gets refactored for thread safety.

Change-Id: Icf39eb4d105555b2eafdfd70b219a73c294900af

src/gromacs/trajectoryanalysis/modules/sasa.cpp
src/gromacs/trajectoryanalysis/modules/surfacearea.cpp
src/gromacs/trajectoryanalysis/modules/surfacearea.h
src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_SurfacePoints12.xml [new file with mode: 0644]
src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_SurfacePoints122.xml [new file with mode: 0644]
src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_SurfacePoints32.xml [new file with mode: 0644]
src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_SurfacePoints42.xml [new file with mode: 0644]
src/gromacs/trajectoryanalysis/tests/surfacearea.cpp

index 32665424fca66e3c30c7ad7b58a6851c2084a662..954ec8f79ac36d9e8d8599851787b295d1b30724 100644 (file)
@@ -387,6 +387,8 @@ class Sasa : public TrajectoryAnalysisModule
          * Empty if the free energy output has not been requested.
          */
         std::vector<real>       dgsFactor_;
+        //! Calculation algorithm.
+        SurfaceAreaCalculator   calculator_;
 
         // Copy and assign disallowed by base.
 };
@@ -609,6 +611,8 @@ Sasa::initAnalysis(const TrajectoryAnalysisSettings &settings,
         }
     }
 
+    calculator_.setDotCount(ndots_);
+
     // Initialize all the output data objects and initialize the output plotters.
 
     area_.setColumnCount(0, 1 + outputSel_.size());
@@ -907,12 +911,10 @@ Sasa::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
     real  totarea, totvolume;
     real *area = NULL, *surfacedots = NULL;
     int   nsurfacedots;
-    int   retval = nsc_dclm_pbc(surfaceSel.coordinates().data(), &radii_[0],
-                                frameData.index_.size(), ndots_, flag, &totarea,
-                                &area, &totvolume, &surfacedots, &nsurfacedots,
-                                &frameData.index_[0],
-                                pbc != NULL ? pbc->ePBC : epbcNONE,
-                                pbc != NULL ? pbc->box : NULL);
+    calculator_.calculate(surfaceSel.coordinates().data(), &radii_[0], pbc,
+                          frameData.index_.size(), &frameData.index_[0], flag,
+                          &totarea, &totvolume, &area,
+                          &surfacedots, &nsurfacedots);
     // Unpack the atomwise areas into the frameData.atomAreas_ array for easier
     // indexing in the case of dynamic surfaceSel.
     if (area != NULL)
@@ -934,10 +936,6 @@ Sasa::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
         sfree(area);
     }
     scoped_guard_sfree dotsGuard(surfacedots);
-    if (retval != 0)
-    {
-        GMX_THROW(InternalError("nsc_dclm_pbc failed"));
-    }
 
     if (bConnolly)
     {
index 9b00a1666816a3ea99f64ffd4f0220136cbdcc96..5f333230e12e6a2e739529707699299ee5291972 100644 (file)
 #include "gromacs/legacyheaders/macros.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/pbcutil/pbc.h"
+#include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/smalloc.h"
 
-#define TEST_NSC 0
-
-#define TEST_ARC 0
-#define TEST_DOD 0
-#define TEST_CUBE 0
-
 #define UNSP_ICO_DOD      9
 #define UNSP_ICO_ARC     10
 
-real   *xpunsp = NULL;
-real    del_cube;
-int     n_dot, ico_cube, last_n_dot = 0, last_densit = 0, last_unsp = 0;
-int     last_cubus = 0;
+static real   *xpunsp = NULL;
+static real    del_cube;
+static int     n_dot, ico_cube, last_n_dot = 0, last_densit = 0, last_unsp = 0;
+static int     last_cubus = 0;
 
 #define FOURPI (4.*M_PI)
 #define TORAD(A)     ((A)*0.017453293)
 #define DP_TOL     0.001
 
 #define UPDATE_FL  __file__ = __FILE__, __line__ = __LINE__
-const char * __file__;   /* declared versions of macros */
-int          __line__;   /* __FILE__  and __LINE__ */
+static const char * __file__;   /* declared versions of macros */
+static int          __line__;   /* __FILE__  and __LINE__ */
 
 #ifdef ERROR
 #undef ERROR
 #endif
 #define ERROR UPDATE_FL, error
-void error(const char *fmt, ...)
+static void error(const char *fmt, ...)
 {
     va_list args;
     fprintf(stderr,
@@ -89,23 +84,8 @@ void error(const char *fmt, ...)
     fprintf(stderr, "\n---> Execution stopped !\n\n");
 }
 
-#define WARNING UPDATE_FL, warning2
-void warning2(const char *fmt, ...)
-{
-    va_list args;
-    fprintf(stderr,
-            "\n---> WARNING : line %i in file %s\n",
-            __line__, __file__);
-    va_start(args, fmt);
-    vfprintf(stderr, fmt, args);
-    va_end(args);
-    fprintf(stderr, " ...!\n\n");
-    fflush(stderr);
-    fflush(stdout);
-}
-
 #define ASIN safe_asin
-real safe_asin(real f)
+static real safe_asin(real f)
 {
     if ( (fabs(f) < 1.00) )
     {
@@ -118,16 +98,13 @@ real safe_asin(real f)
     return(M_PI_2);
 }
 
-
-
-
 /* routines for dot distributions on the surface of the unit sphere */
-real rg, rh;
+static real rh;
 
-void icosaeder_vertices(real *xus)
+static void icosaeder_vertices(real *xus)
 {
     rh = sqrt(1.-2.*cos(TORAD(72.)))/(1.-cos(TORAD(72.)));
-    rg = cos(TORAD(72.))/(1.-cos(TORAD(72.)));
+    const real rg = cos(TORAD(72.))/(1.-cos(TORAD(72.)));
     /* icosaeder vertices */
     xus[ 0] = 0.;                  xus[ 1] = 0.;                  xus[ 2] = 1.;
     xus[ 3] = rh*cos(TORAD(72.));  xus[ 4] = rh*sin(TORAD(72.));  xus[ 5] = rg;
@@ -144,9 +121,9 @@ void icosaeder_vertices(real *xus)
 }
 
 
-void divarc(real x1, real y1, real z1,
-            real x2, real y2, real z2,
-            int div1, int div2, real *xr, real *yr, real *zr)
+static void divarc(real x1, real y1, real z1,
+                   real x2, real y2, real z2,
+                   int div1, int div2, real *xr, real *yr, real *zr)
 {
 
     real xd, yd, zd, dd, d1, d2, s, x, y, z;
@@ -184,7 +161,7 @@ void divarc(real x1, real y1, real z1,
     *xr = x/dd; *yr = y/dd; *zr = z/dd;
 }
 
-int ico_dot_arc(int densit)   /* densit...required dots per unit sphere */
+static int ico_dot_arc(int densit)   /* densit...required dots per unit sphere */
 {
     /* dot distribution on a unit sphere based on an icosaeder *
      * great circle average refining of icosahedral face       */
@@ -318,9 +295,9 @@ int ico_dot_arc(int densit)   /* densit...required dots per unit sphere */
     }   /* end of if (tess > 1) */
 
     return n_dot;
-}                           /* end of routine ico_dot_arc */
+}                                  /* end of routine ico_dot_arc */
 
-int ico_dot_dod(int densit) /* densit...required dots per unit sphere */
+static int ico_dot_dod(int densit) /* densit...required dots per unit sphere */
 {
     /* dot distribution on a unit sphere based on an icosaeder *
      * great circle average refining of icosahedral face       */
@@ -505,7 +482,7 @@ int ico_dot_dod(int densit) /* densit...required dots per unit sphere */
     return n_dot;
 }       /* end of routine ico_dot_dod */
 
-int unsp_type(int densit)
+static int unsp_type(int densit)
 {
     int i1, i2;
     i1 = 1;
@@ -528,10 +505,10 @@ int unsp_type(int densit)
     }
 }
 
-int make_unsp(int densit, int mode, int * num_dot, int cubus)
+static void make_unsp(int densit, int mode, int * num_dot, int cubus)
 {
     int  *ico_wk, *ico_pt;
-    int   ndot, ico_cube_cb, i, j, k, l, ijk, tn, tl, tl2;
+    int   ico_cube_cb, i, j, k, l, ijk, tn, tl, tl2;
     real *xus;
     int  *work;
     real  x, y, z;
@@ -541,10 +518,14 @@ int make_unsp(int densit, int mode, int * num_dot, int cubus)
         free(xpunsp);
     }
 
-    k = 1; if (mode < 0)
+    k = 1;
+    if (mode < 0)
     {
-        k = 0; mode = -mode;
+        k    = 0;
+        mode = -mode;
     }
+    // Assignment to satisfy cppcheck.
+    int ndot = 0;
     if (mode == UNSP_ICO_ARC)
     {
         ndot = ico_dot_arc(densit);
@@ -555,14 +536,14 @@ int make_unsp(int densit, int mode, int * num_dot, int cubus)
     }
     else
     {
-        WARNING("make_unsp: mode %c%d not allowed", (k) ? '+' : '-', mode);
-        return 1;
+        GMX_RELEASE_ASSERT(false, "Invalid unit sphere mode");
     }
 
     last_n_dot = ndot; last_densit = densit; last_unsp = mode;
-    *num_dot   = ndot; if (k)
+    *num_dot   = ndot;
+    if (k)
     {
-        return 0;
+        return;
     }
 
     /* in the following the dots of the unit sphere may be resorted */
@@ -646,8 +627,6 @@ int make_unsp(int densit, int mode, int * num_dot, int cubus)
     }         /* cycle i */
     sfree(ico_wk);
     sfree(work);
-
-    return 0;
 }
 
 
@@ -658,12 +637,14 @@ typedef struct _stwknb {
     real dot;
 } Neighb;
 
-int nsc_dclm_pbc(const rvec *coords, real *radius, int nat,
-                 int  densit, int mode,
-                 real *value_of_area, real **at_area,
-                 real *value_of_vol,
-                 real **lidots, int *nu_dots,
-                 atom_id index[], int ePBC, matrix box)
+// This is extern "C" for now to make it easy to put valgrind suppressions.
+extern "C" void
+nsc_dclm_pbc(const rvec *coords, real *radius, int nat,
+             int  densit, int mode,
+             real *value_of_area, real **at_area,
+             real *value_of_vol,
+             real **lidots, int *nu_dots,
+             atom_id index[], int ePBC, matrix box)
 {
     int         iat, i, ii, iii, ix, iy, iz, ixe, ixs, iye, iys, ize, izs, i_ac;
     int         jat, j, jj, jjj, jx, jy, jz;
@@ -689,10 +670,7 @@ int nsc_dclm_pbc(const rvec *coords, real *radius, int nat,
     if (distribution != -last_unsp || last_cubus != 4 ||
         (densit != last_densit && densit != last_n_dot))
     {
-        if (make_unsp(densit, (-distribution), &n_dot, 4))
-        {
-            return 1;
-        }
+        make_unsp(densit, (-distribution), &n_dot, 4);
     }
     xus = xpunsp;
 
@@ -708,8 +686,7 @@ int nsc_dclm_pbc(const rvec *coords, real *radius, int nat,
     /* calculate neighbour list with the box algorithm */
     if (nat == 0)
     {
-        WARNING("nsc_dclm: no surface atoms selected");
-        return 1;
+        return;
     }
     if (mode & FLAG_VOLUME)
     {
@@ -1144,86 +1121,115 @@ int nsc_dclm_pbc(const rvec *coords, real *radius, int nat,
     {
         fprintf(debug, "area=%8.3f\n", area);
     }
-
-    return 0;
 }
 
+namespace gmx
+{
+
+class SurfaceAreaCalculator::Impl
+{
+    public:
+        Impl() : dotCount_(32), flags_(0)
+        {
+        }
+
+        int dotCount_;
+        int flags_;
+};
+
+SurfaceAreaCalculator::SurfaceAreaCalculator()
+    : impl_(new Impl())
+{
+}
 
-#if TEST_NSC > 0
-#define NAT 2
-main () {
+SurfaceAreaCalculator::~SurfaceAreaCalculator()
+{
+}
 
-    int    i, j, ndots;
-    real   co[3*NAT], ra[NAT], area, volume, a, b, c;
-    real * dots;
-    real * at_area;
-    FILE  *fp;
+void SurfaceAreaCalculator::setDotCount(int dotCount)
+{
+    impl_->dotCount_ = dotCount;
+}
 
+void SurfaceAreaCalculator::setCalculateVolume(bool bVolume)
+{
+    if (bVolume)
+    {
+        impl_->flags_ |= FLAG_VOLUME;
+    }
+    else
+    {
+        impl_->flags_ &= ~FLAG_VOLUME;
+    }
+}
 
-    a  = 1.; c = 0.1;
-    fp = fopen("nsc.txt", "w+");
-    for (i = 1; i <= NAT; i++)
+void SurfaceAreaCalculator::setCalculateAtomArea(bool bAtomArea)
+{
+    if (bAtomArea)
     {
-        j         = i-1;
-        co[3*i-3] = j*1*c;
-        co[3*i-2] = j*1*c;
-        co[3*i-1] = j*1*c;
-        /*
-           co[3*i-3] = i*1.4;
-           co[3*i-2] = 0.;
-           co[3*i-1] = 0.;
-         */
-        /*
-           co[3*i-2] = a*0.3;
-           a = -a; b=0;
-           if (i%3 == 0) b=0.5;
-           co[3*i-1] = b;
-           ra[i-1] = 2.0;
-         */
-        ra[i-1] = (1.+j*0.5)*c;
+        impl_->flags_ |= FLAG_ATOM_AREA;
     }
-    /*
-       if (NSC(co, ra, NAT, 42, NULL, &area,
-     */
-    if (NSC(co, ra, NAT, 42, NULL, &area,
-            NULL, NULL, NULL, NULL))
+    else
     {
-        ERROR("error in NSC");
+        impl_->flags_ &= ~FLAG_ATOM_AREA;
     }
-    fprintf(fp, "\n");
-    fprintf(fp, "area     : %8.3f\n", area);
-    fprintf(fp, "\n");
-    fprintf(fp, "\n");
-    fprintf(fp, "next call\n");
-    fprintf(fp, "\n");
-    fprintf(fp, "\n");
-
-    if (NSC(co, ra, NAT, 42, FLAG_VOLUME | FLAG_ATOM_AREA | FLAG_DOTS, &area,
-            &at_area, &volume,
-            &dots, &ndots))
+}
+
+void SurfaceAreaCalculator::setCalculateSurfaceDots(bool bDots)
+{
+    if (bDots)
+    {
+        impl_->flags_ |= FLAG_DOTS;
+    }
+    else
     {
-        ERROR("error in NSC");
+        impl_->flags_ &= ~FLAG_DOTS;
     }
+}
 
-    fprintf(fp, "\n");
-    fprintf(fp, "area     : %8.3f\n", area);
-    printf("area     : %8.3f\n", area);
-    fprintf(fp, "volume   : %8.3f\n", volume);
-    printf("volume   : %8.3f\n", volume);
-    fprintf(fp, "ndots    : %8d\n", ndots);
-    printf("ndots    : %8d\n", ndots);
-    fprintf(fp, "\n");
-    for (i = 1; i <= NAT; i++)
+void SurfaceAreaCalculator::calculate(
+        const rvec *x, const real *radius, const t_pbc *pbc,
+        int nat, atom_id index[], int flags, real *area, real *volume,
+        real **at_area, real **lidots, int *n_dots) const
+{
+    flags |= impl_->flags_;
+    *area  = 0;
+    if (volume == NULL)
     {
-        fprintf(fp, "%4d ATOM %7.2f %7.2f %7.2f  ra=%4.1f  area=%8.3f\n",
-                i, co[3*i-3], co[3*i-2], co[3*i-1], ra[i-1], at_area[i-1]);
+        flags &= ~FLAG_VOLUME;
     }
-    fprintf(fp, "\n");
-    fprintf(fp, "DOTS : %8d\n", ndots);
-    for (i = 1; i <= ndots; i++)
+    else
     {
-        fprintf(fp, "%4d DOTS %8.2f %8.2f %8.2f\n",
-                i, dots[3*i-3], dots[3*i-2], dots[3*i-1]);
+        *volume = 0;
     }
+    if (at_area == NULL)
+    {
+        flags &= ~FLAG_ATOM_AREA;
+    }
+    else
+    {
+        *at_area = NULL;
+    }
+    if (lidots == NULL)
+    {
+        flags &= ~FLAG_DOTS;
+    }
+    else
+    {
+        *lidots = NULL;
+    }
+    if (n_dots == NULL)
+    {
+        flags &= ~FLAG_DOTS;
+    }
+    else
+    {
+        *n_dots = 0;
+    }
+    nsc_dclm_pbc(x, const_cast<real *>(radius), nat, impl_->dotCount_,
+                 flags, area, at_area, volume, lidots, n_dots, index,
+                 pbc != NULL ? pbc->ePBC : epbcNONE,
+                 pbc != NULL ? const_cast<rvec *>(pbc->box) : NULL);
 }
-#endif
+
+} // namespace gmx
index eb5129ea22d8de6d0da3ea67c35a676a78a02371..25270eda8d9f8caf458bd5fc9549097e7f6338dc 100644 (file)
 #define GMX_TRAJECTORYANALYSIS_SURFACEAREA_H
 
 #include "gromacs/legacyheaders/types/simple.h"
+#include "gromacs/utility/classhelpers.h"
+
+struct t_pbc;
 
 #define FLAG_DOTS       01
 #define FLAG_VOLUME     02
 #define FLAG_ATOM_AREA  04
 
-#ifdef __cplusplus
-extern "C"
+namespace gmx
 {
-#endif
-
-int nsc_dclm_pbc(const rvec *coords, real *radius, int nat,
-                 int  densit, int mode,
-                 real *value_of_area, real **at_area,
-                 real *value_of_vol,
-                 real **lidots, int *nu_dots,
-                 atom_id index[], int ePBC, matrix box);
-
-/*
-    User notes :
-   The input requirements :
-   The arrays with atom coordinates and radii are thought to start
-   with index 0, i.e., places 0, 1, and 2 are the x-, y-, and z-
-   coordinates of the zero-th atom and place 0 in the other array
-   is its radius.
-
-   PLEASE TAKE INTO ACCOUNT THAT THE RADII GIVEN HERE ARE DIRECTLY
-   USED FOR SURFACE CALCULATION. NSC does not increment with a probe
-   radius.
-
-   The user can define any number of dots. The program selects a
-   dot density that is the lowest possible with at least the required
-   number of dots. The points are distributed in accordance with the
-   icosahedron-based or the dodecahedron-based method as described in
-   ref. 1.
-
-   The output requirements are :
-   1 and 3 :  pointer to an existing real
-   2 and 4 :  pointer to an existing pointer to real
-             NSC allocates memory for an array
-   5       :  pointer to an existing integer
-
-   The subroutine NSC makes use of variant 2 described in reference 1.
-   By selecting the necessary output via flags, the requirements for
-   cpu-time and computer memory can be adapted to the actual needs.
-
-   Example : flag = FLAG_VOLUME | FLAG_ATOM_AREA | FLAG_DOTS
-          The routine calculates the area, volume and the dot surface. The
-          program allocates arrays for the atomwise areas and for the surface
-          dots. The addresses are returned in the pointers to pointers to
-          real.
-          This variant is not recommended because normally the dot surface
-          is needed for low point density (e.g.42) at which area and volume
-          are inaccurate. The sign "|" is used as binary AND !
-
-          flag = FLAG_VOLUME | FLAG_ATOM_AREA
-          In this case the large arrays for storing the surface dots
-          are not allocated. A large point number of the fully accessible
-          sphere can be selected. Good accuracy is already achieved with
-          600-700 points per sphere (accuracy of about 1.5 square Angstrem
-          per atomic sphere).
-          Output pointers 4 and 5 may be NULL.
-
-          flag = FLAG_DOTS
-          Only the dot surface is produced.
-          Output pointers 2 and 3 may be NULL.
-
-   The output pointer 1 cannot be set to NULL in any circumstances. The
-   overall area value is returned in every mode.
-
-   All files calling NSC should include surfacearea.h !!
-
-
-   Example for calling NSC (contents of user file):
-
-   ...
-   #include "surfacearea.h"
-
-   int routine_calling_NSC(int n_atom, real *coordinates, real *radii) {
-   real area, volume, *atomwise_area, *surface_dots;
-   int    i, density = 300, n_dots;
-
-   ...
-
-   for (i=0; i<n_atom; i++) {
-   radii[i]  += 1.4      /# add the probe radius if necessary #/
-
-   if (NSC(coordinates, radii, n_atom, density,
-          FLAG_AREA | FLAG_VOLUME | FLAG_DOTS,
-          &area, &atomwise_area, &volume, &surface_dots, &n_dots))
-    printf("error occured\n");
-    return 1;
-    }
-
-   ...
-
-   /# do something with areas, volume and surface dots #/
-
-   ...
-
-   return 0;
-   }
 
+/*! \internal
+ * \brief
+ * Computes surface areas for a group of atoms/spheres.
+ *
+ * This class provides a surface area/volume calculator.
+ *
+ * The algorithm is based on representing each atom/sphere surface as a set of
+ * dots, and determining which dots are on the surface (not covered by any
+ * other atom/sphere).  The dots are distributed evenly using an icosahedron- or
+ * a dodecahedron-based method (see the original reference cited in the code).
+ * The area is then estimated from the area represented by each dot.
+ * The volume is calculated by selecting a fixed point and integrating over the
+ * surface dots, summing up the cones whose tip is at the fixed point and base
+ * at the surface points.
+ *
+ * The default dot density per sphere is 32, which gives quite inaccurate
+ * areas and volumes, but a reasonable number of surface points.  According to
+ * original documentation of the method, a density of 600-700 dots gives an
+ * accuracy of 1.5 A^2 per atom.
+ *
+ * \ingroup module_trajectoryanalysis
  */
-
-#ifdef __cplusplus
-}
-#endif
+class SurfaceAreaCalculator
+{
+    public:
+        /*! \brief
+         * Initializes a surface area calculator.
+         *
+         * \throws std::bad_alloc if out of memory.
+         */
+        SurfaceAreaCalculator();
+        ~SurfaceAreaCalculator();
+
+        /*! \brief
+         * Sets the number of surface dots per sphere to use.
+         *
+         * If not called, the default is 32.
+         */
+        void setDotCount(int dotCount);
+
+        /*! \brief
+         * Requests calculation of volume.
+         *
+         * If not called, and FLAG_VOLUME is not passed to calculate(), the
+         * volume output is not produced.
+         *
+         * Does not throw.
+         */
+        void setCalculateVolume(bool bVolume);
+        /*! \brief
+         * Requests output of per-atom areas.
+         *
+         * If not called, and FLAG_ATOM_AREA is not passed to calculate(), the
+         * atom area output is not produced.
+         *
+         * Does not throw.
+         */
+        void setCalculateAtomArea(bool bAtomArea);
+        /*! \brief
+         * Requests output of all surface dots.
+         *
+         * If not called, and FLAG_DOTS is not passed to calculate(), the
+         * surface dot output is not produced.
+         *
+         * Does not throw.
+         */
+        void setCalculateSurfaceDots(bool bDots);
+
+        /*! \brief
+         * Calculates the surface area for a set of positions.
+         *
+         * \param[in]  x       Atom positions (sphere centers).
+         * \param[in]  radius  Radius for each atom/sphere.
+         *     These radii are used as-is, without adding any probe radius.
+         * \param[in]  pbc     PBC information (if `NULL`, calculation is done
+         *     without PBC).
+         * \param[in]  nat     Number of atoms to calculate.
+         * \param[in]  index   Atom indices to include in the calculation.
+         * \param[in]  flags   Additional flags for the calculation.
+         * \param[out] area    Total surface area (must be non-`NULL`).
+         * \param[out] volume  Total volume (can be `NULL`).
+         * \param[out] at_area Surface area for each atom in \p index
+         *     (\p nat values) (can be `NULL`).
+         * \param[out] lidots  Surface dots as x,y,z triplets (`3*lidots` values)
+         *     (can be `NULL`).
+         * \param[out] n_dots Number of surface dots in \p lidots
+         *     (can be `NULL`).
+         *
+         * Calculates the surface area of spheres centered at `x[index[0]]`,
+         * ..., `x[index[nat-1]]`, with radii `radius[index[0]]`, ....
+         *
+         * If \p flags is 0, the calculation is done for the items specified
+         * with setCalculateVolume(), setCalculateAtomArea(), and
+         * setCalculateSurfaceDots().  Flags can specify FLAG_VOLUME,
+         * FLAG_ATOM_AREA, and/or FLAG_DOTS to request additional output for
+         * this particular calculation.  If any output is `NULL`, that output
+         * is not calculated, irrespective of the calculation mode set.
+         *
+         * \todo
+         * Make the output options more C++-like, in particular for the array
+         * outputs.
+         *
+         * \todo
+         * Make this thread-safe; currently, only a single concurrent call to
+         * calculate() is supported.
+         */
+        void calculate(const rvec *x, const real *radius, const t_pbc *pbc,
+                       int nat, atom_id index[], int flags, real *area,
+                       real *volume, real **at_area,
+                       real **lidots, int *n_dots) const;
+
+    private:
+        class Impl;
+
+        PrivateImplPointer<Impl> impl_;
+};
+
+} // namespace gmx
 
 #endif
diff --git a/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_SurfacePoints12.xml b/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_SurfacePoints12.xml
new file mode 100644 (file)
index 0000000..bf3932b
--- /dev/null
@@ -0,0 +1,46 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <SASA Name="Surface">
+    <Real Name="Area">12.566370614359172</Real>
+    <Sequence Name="Dots">
+      <Int Name="Length">36</Int>
+      <Real>0</Real>
+      <Real>0</Real>
+      <Real>-1</Real>
+      <Real>-0.27639310998086691</Real>
+      <Real>-0.85065087452807109</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>0.72360690737127242</Real>
+      <Real>-0.52573101980466264</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>-0.8944272254243314</Real>
+      <Real>0</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>0.72360681651420888</Real>
+      <Real>0.52573114485870043</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>-0.27639325699069922</Real>
+      <Real>0.85065082676168668</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>0.27639333049561254</Real>
+      <Real>-0.85065080287848482</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>-0.72360677108566906</Real>
+      <Real>-0.52573120738571333</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>0.8944272254243314</Real>
+      <Real>0</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>-0.72360686194274337</Real>
+      <Real>0.52573108233168353</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>0.27639318348578412</Real>
+      <Real>0.85065085064488211</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>0</Real>
+      <Real>0</Real>
+      <Real>1</Real>
+    </Sequence>
+  </SASA>
+</ReferenceData>
diff --git a/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_SurfacePoints122.xml b/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_SurfacePoints122.xml
new file mode 100644 (file)
index 0000000..71fd5e3
--- /dev/null
@@ -0,0 +1,376 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <SASA Name="Surface">
+    <Real Name="Area">12.566370614359172</Real>
+    <Sequence Name="Dots">
+      <Int Name="Length">366</Int>
+      <Real>0</Real>
+      <Real>0</Real>
+      <Real>-1</Real>
+      <Real>0.099017100750074208</Real>
+      <Real>-0.30474315751655773</Real>
+      <Real>-0.94727357278976609</Real>
+      <Real>-0.25923000773171251</Real>
+      <Real>-0.18834164910928383</Real>
+      <Real>-0.94727357521584643</Real>
+      <Real>0.32042594696241011</Real>
+      <Real>2.7688072479927719e-08</Real>
+      <Real>-0.94727356793760586</Real>
+      <Real>-0.25923003561123026</Real>
+      <Real>0.18834163514060129</Real>
+      <Real>-0.94727357036368598</Real>
+      <Real>0.099017048084233991</Real>
+      <Real>0.30474317462872957</Real>
+      <Real>-0.94727357278976609</Real>
+      <Real>-0.16245980125708112</Real>
+      <Real>-0.50000003709380847</Real>
+      <Real>-0.85065079549818778</Real>
+      <Real>0.42532548914805512</Real>
+      <Real>-0.3090169737533614</Real>
+      <Real>-0.85065077335724681</Real>
+      <Real>-0.52573114485869987</Real>
+      <Real>-1.5142845699589461e-08</Real>
+      <Real>-0.85065078811787431</Real>
+      <Real>0.42532543574360737</Real>
+      <Real>0.30901704725828416</Real>
+      <Real>-0.8506507733572467</Real>
+      <Real>-0.1624598806459048</Real>
+      <Real>0.50000003641109625</Real>
+      <Real>-0.85065078073756029</Real>
+      <Real>0.18759256558961399</Real>
+      <Real>-0.57735027920788806</Real>
+      <Real>-0.79465444341177649</Real>
+      <Real>-0.49112347245450155</Real>
+      <Real>-0.35682213462759943</Real>
+      <Real>-0.79465445260442358</Real>
+      <Real>0.60706206007773655</Real>
+      <Real>5.2456358414752335e-08</Real>
+      <Real>-0.79465442502648409</Real>
+      <Real>-0.49112352275791088</Real>
+      <Real>0.35682210633554395</Real>
+      <Real>-0.79465443421912996</Real>
+      <Real>0.18759246581169672</Real>
+      <Real>0.57735031162770434</Real>
+      <Real>-0.79465444341177638</Real>
+      <Real>-0.046871646661550442</Real>
+      <Real>-0.75374273850500617</Real>
+      <Real>-0.65549594422102153</Real>
+      <Real>-0.40511875788237067</Real>
+      <Real>-0.6373412338447112</Real>
+      <Real>-0.6554959524307884</Real>
+      <Real>0.48095898541609272</Real>
+      <Real>-0.58224008918771142</Real>
+      <Real>-0.65549594422102164</Real>
+      <Real>0.70236782279174126</Real>
+      <Real>-0.2774969011402712</Real>
+      <Real>-0.65549592780148958</Real>
+      <Real>-0.73133609169283587</Real>
+      <Real>-0.18834165055637986</Real>
+      <Real>-0.6554959524307884</Real>
+      <Real>-0.73133611075236105</Real>
+      <Real>0.18834163369350529</Real>
+      <Real>-0.65549593601125533</Real>
+      <Real>0.70236777483460133</Real>
+      <Real>0.27749702252376535</Real>
+      <Real>-0.65549592780148958</Real>
+      <Real>0.48095888479311216</Real>
+      <Real>0.58224017230723379</Real>
+      <Real>-0.65549594422102153</Real>
+      <Real>-0.40511885788011226</Real>
+      <Real>0.63734118716951105</Real>
+      <Real>-0.65549593601125544</Real>
+      <Real>-0.046871776923682162</Real>
+      <Real>0.75374273040461759</Real>
+      <Real>-0.65549594422102153</Real>
+      <Real>0.26286567453986009</Real>
+      <Real>-0.80901697843275555</Real>
+      <Real>-0.52573107741148295</Real>
+      <Real>-0.68819093826563327</Real>
+      <Real>-0.50000005836206629</Real>
+      <Real>-0.52573108537254831</Real>
+      <Real>0.85065083964296051</Real>
+      <Real>5.6152759798857928e-08</Real>
+      <Real>-0.52573106148935433</Real>
+      <Real>-0.68819099963114494</Real>
+      <Real>0.49999999064120482</Real>
+      <Real>-0.52573106945041825</Real>
+      <Real>0.26286553472520674</Real>
+      <Real>0.80901702386129837</Real>
+      <Real>-0.52573107741148284</Real>
+      <Real>-0.27639310998086691</Real>
+      <Real>-0.85065087452807109</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>0.72360690737127242</Real>
+      <Real>-0.52573101980466264</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>-0.8944272254243314</Real>
+      <Real>0</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>0.72360681651420888</Real>
+      <Real>0.52573114485870043</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>-0.27639325699069922</Real>
+      <Real>0.85065082676168668</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>0.014324278037483651</Real>
+      <Real>-0.9420843338704582</Real>
+      <Real>-0.33507002691162935</Real>
+      <Real>0.54215490657597232</Real>
+      <Real>-0.77058168570310492</Real>
+      <Real>-0.33507002691162935</Real>
+      <Real>-0.56533168971421621</Real>
+      <Real>-0.75374276489632219</Real>
+      <Real>-0.33507002995092922</Real>
+      <Real>-0.89154902505459532</Real>
+      <Real>-0.30474318677637491</Real>
+      <Real>-0.33507003154203469</Real>
+      <Real>0.90040188534718624</Real>
+      <Real>-0.27749689752267748</Real>
+      <Real>-0.33507001765082056</Real>
+      <Real>0.90040184594129713</Real>
+      <Real>0.27749702154159639</Real>
+      <Real>-0.33507002083303039</Real>
+      <Real>-0.89154904724027939</Real>
+      <Real>0.30474313205289755</Real>
+      <Real>-0.33507002228122468</Real>
+      <Real>-0.56533180027257979</Real>
+      <Real>0.75374268467611705</Real>
+      <Real>-0.33507002387232976</Real>
+      <Real>0.54215477340371809</Real>
+      <Real>0.77058177939852901</Real>
+      <Real>-0.33507002691162935</Real>
+      <Real>0.014324115226078571</Real>
+      <Real>0.94208433634597177</Real>
+      <Real>-0.33507002691162935</Real>
+      <Real>0.30353113023123213</Real>
+      <Real>-0.93417232306880937</Real>
+      <Real>-0.1875924406599985</Real>
+      <Real>-0.79465443246348733</Real>
+      <Real>-0.57735033372812661</Real>
+      <Real>-0.18759244417127635</Real>
+      <Real>0.98224695410165486</Real>
+      <Real>5.2456354487327486e-08</Real>
+      <Real>-0.18759243363744385</Real>
+      <Real>-0.79465450249337599</Real>
+      <Real>0.57735023962202381</Real>
+      <Real>-0.18759243714872095</Real>
+      <Real>0.30353096878717961</Real>
+      <Real>0.93417237552517096</Real>
+      <Real>-0.18759244065999844</Real>
+      <Real>-0.30610163692434295</Real>
+      <Real>-0.94208435645226973</Real>
+      <Real>-0.13703595586615472</Real>
+      <Real>0.80138486174169055</Real>
+      <Real>-0.58224002713709733</Real>
+      <Real>-0.13703595940717078</Real>
+      <Real>-0.99056607391927598</Real>
+      <Real>-2.7688068789135795e-08</Real>
+      <Real>-0.13703595586615464</Real>
+      <Real>0.80138477778966677</Real>
+      <Real>0.58224014352057152</Real>
+      <Real>-0.13703595586615475</Real>
+      <Real>-0.30610179973574703</Real>
+      <Real>0.94208430355164396</Real>
+      <Real>-0.13703595586615469</Real>
+      <Real>1.2961531081767674e-07</Real>
+      <Real>-0.99999999999999167</Real>
+      <Real>-5.5511151231257827e-17</Real>
+      <Real>-0.58778517212449111</Real>
+      <Real>-0.80901705262038981</Real>
+      <Real>-1.8793522849058775e-09</Real>
+      <Real>0.58778534971183316</Real>
+      <Real>-0.80901692359563038</Real>
+      <Real>-3.7587038759223645e-09</Real>
+      <Real>-0.95105649618237209</Real>
+      <Real>-0.30901705627571729</Real>
+      <Real>-1.879352201639151e-09</Real>
+      <Real>0.95105654032715115</Real>
+      <Real>-0.30901692041205481</Real>
+      <Real>2.7755575615628914e-17</Real>
+      <Real>-0.95105653106748755</Real>
+      <Real>0.3090169489103749</Real>
+      <Real>1.879352062861273e-09</Real>
+      <Real>0.95105650544203846</Real>
+      <Real>0.30901702777739931</Real>
+      <Real>3.7587038759223654e-09</Real>
+      <Real>-0.58778529051605588</Real>
+      <Real>0.80901696660388756</Real>
+      <Real>1.8793520073501213e-09</Real>
+      <Real>0.58778523132027516</Real>
+      <Real>0.80901700961214085</Real>
+      <Real>2.7755575615628914e-17</Real>
+      <Real>-4.3205103550381097e-08</Real>
+      <Real>0.99999999999999911</Real>
+      <Real>-2.7755575615628914e-17</Real>
+      <Real>0.30610188114144588</Real>
+      <Real>-0.94208427710132048</Real>
+      <Real>0.13703595586615475</Real>
+      <Real>-0.80138473581364844</Real>
+      <Real>-0.58224020171230428</Real>
+      <Real>0.13703595409564634</Real>
+      <Real>0.99056607342940817</Real>
+      <Real>2.7688069489002942e-08</Real>
+      <Real>0.13703595940717081</Real>
+      <Real>-0.80138481976568088</Real>
+      <Real>0.58224008532883575</Real>
+      <Real>0.13703595763666276</Real>
+      <Real>0.30610171833004618</Real>
+      <Real>0.94208433000196035</Real>
+      <Real>0.13703595586615469</Real>
+      <Real>-0.30353088806515</Real>
+      <Real>-0.93417240175334149</Real>
+      <Real>0.1875924406599985</Real>
+      <Real>0.79465453750831549</Real>
+      <Real>-0.5773501925689698</Real>
+      <Real>0.18759243363744388</Real>
+      <Real>-0.98224695276046647</Real>
+      <Real>-5.245635309796613e-08</Real>
+      <Real>0.18759244065999858</Real>
+      <Real>0.79465446747843327</Real>
+      <Real>0.57735028667507615</Real>
+      <Real>0.18759244065999844</Real>
+      <Real>-0.3035310495092069</Real>
+      <Real>0.93417234929699366</Real>
+      <Real>0.18759244065999844</Real>
+      <Real>-0.014324033820375781</Real>
+      <Real>-0.94208433758371801</Real>
+      <Real>0.33507002691162935</Real>
+      <Real>-0.54215470681758493</Real>
+      <Real>-0.77058182624623239</Real>
+      <Real>0.33507002691162935</Real>
+      <Real>0.56533185555175702</Real>
+      <Real>-0.75374264456600837</Real>
+      <Real>0.33507002083303045</Real>
+      <Real>0.8915490583331197</Real>
+      <Real>-0.30474310469115867</Real>
+      <Real>0.33507001765082062</Real>
+      <Real>-0.90040183190415546</Real>
+      <Real>-0.27749705782731754</Real>
+      <Real>0.33507002850273465</Real>
+      <Real>-0.90040187131004945</Real>
+      <Real>0.27749693380840212</Real>
+      <Real>0.33507002532052427</Real>
+      <Real>0.89154903614743786</Real>
+      <Real>0.30474315941463637</Real>
+      <Real>0.33507002691162935</Real>
+      <Real>0.56533174499339967</Real>
+      <Real>0.75374272478622162</Real>
+      <Real>0.33507002691162935</Real>
+      <Real>-0.54215483998984715</Real>
+      <Real>0.77058173255081996</Real>
+      <Real>0.33507002691162935</Real>
+      <Real>-0.014324196631781194</Real>
+      <Real>0.94208433510821854</Real>
+      <Real>0.33507002691162929</Real>
+      <Real>0.27639333049561254</Real>
+      <Real>-0.85065080287848482</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>-0.72360677108566906</Real>
+      <Real>-0.52573120738571333</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>0.8944272254243314</Real>
+      <Real>0</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>-0.72360686194274337</Real>
+      <Real>0.52573108233168353</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>0.27639318348578412</Real>
+      <Real>0.85065085064488211</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>-0.26286546481787715</Real>
+      <Real>-0.80901704657556062</Real>
+      <Real>0.52573107741148295</Real>
+      <Real>0.68819103031389739</Real>
+      <Real>-0.49999995678077147</Real>
+      <Real>0.52573106148935422</Real>
+      <Real>-0.85065082980254492</Real>
+      <Real>-5.6152758158248256e-08</Real>
+      <Real>0.52573107741148306</Real>
+      <Real>0.68819096894839027</Real>
+      <Real>0.5000000245016365</Real>
+      <Real>0.52573107741148284</Real>
+      <Real>-0.26286560463253433</Real>
+      <Real>0.80901700114703001</Real>
+      <Real>0.52573107741148306</Real>
+      <Real>0.046871842054747564</Real>
+      <Real>-0.75374272635441475</Real>
+      <Real>0.65549594422102153</Real>
+      <Real>0.40511890787898003</Real>
+      <Real>-0.63734116383190476</Real>
+      <Real>0.65549592780148969</Real>
+      <Real>-0.48095883448161647</Real>
+      <Real>-0.58224021386698843</Real>
+      <Real>0.65549594422102153</Real>
+      <Real>-0.70236774534829738</Real>
+      <Real>-0.27749705837024763</Real>
+      <Real>0.65549594422102153</Real>
+      <Real>0.7313361202821228</Real>
+      <Real>-0.18834162526206655</Real>
+      <Real>0.65549592780148969</Real>
+      <Real>0.73133610122259873</Real>
+      <Real>0.18834164212494314</Real>
+      <Real>0.65549594422102153</Real>
+      <Real>-0.70236779330544363</Real>
+      <Real>0.27749693698675881</Real>
+      <Real>0.65549594422102153</Real>
+      <Real>-0.48095893510460419</Real>
+      <Real>0.58224013074747483</Real>
+      <Real>0.65549594422102153</Real>
+      <Real>0.4051188078812426</Real>
+      <Real>0.63734121050711312</Real>
+      <Real>0.65549594422102153</Real>
+      <Real>0.046871711792616434</Real>
+      <Real>0.75374273445481466</Real>
+      <Real>0.65549594422102175</Real>
+      <Real>-0.18759241592273596</Real>
+      <Real>-0.57735032783760598</Real>
+      <Real>0.79465444341177638</Real>
+      <Real>0.49112354790961382</Real>
+      <Real>-0.35682209218951294</Real>
+      <Real>0.79465442502648398</Real>
+      <Real>-0.60706203601107944</Real>
+      <Real>-5.2456356275319567e-08</Real>
+      <Real>0.79465444341177627</Real>
+      <Real>0.49112349760620688</Real>
+      <Real>0.3568221204815728</Real>
+      <Real>0.79465444341177638</Real>
+      <Real>-0.18759251570065602</Real>
+      <Real>0.57735029541779825</Real>
+      <Real>0.79465444341177627</Real>
+      <Real>0.16245992034031653</Real>
+      <Real>-0.50000003606973475</Real>
+      <Real>0.8506507733572467</Real>
+      <Real>-0.42532539862021479</Real>
+      <Real>-0.30901705772173071</Real>
+      <Real>0.85065078811787398</Real>
+      <Real>0.52573116874189396</Real>
+      <Real>1.5142847809013205e-08</Real>
+      <Real>0.85065077335724681</Real>
+      <Real>-0.42532545202466426</Real>
+      <Real>0.30901698421681456</Real>
+      <Real>0.85065078811787398</Real>
+      <Real>0.16245984095149318</Real>
+      <Real>0.50000003675245419</Real>
+      <Real>0.85065078811787398</Real>
+      <Real>-0.099017021751312745</Real>
+      <Real>-0.30474318318481203</Real>
+      <Real>0.94727357278976609</Real>
+      <Real>0.25923004955098822</Real>
+      <Real>-0.18834162815625813</Real>
+      <Real>0.94727356793760586</Real>
+      <Real>-0.32042593261799379</Real>
+      <Real>-2.7688071208844714e-08</Real>
+      <Real>0.94727357278976598</Real>
+      <Real>0.25923002167147174</Real>
+      <Real>0.1883416421249432</Real>
+      <Real>0.94727357278976609</Real>
+      <Real>-0.09901707441715446</Real>
+      <Real>0.30474316607264479</Real>
+      <Real>0.94727357278976609</Real>
+      <Real>0</Real>
+      <Real>0</Real>
+      <Real>1</Real>
+    </Sequence>
+  </SASA>
+</ReferenceData>
diff --git a/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_SurfacePoints32.xml b/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_SurfacePoints32.xml
new file mode 100644 (file)
index 0000000..20dab27
--- /dev/null
@@ -0,0 +1,106 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <SASA Name="Surface">
+    <Real Name="Area">12.566370614359172</Real>
+    <Sequence Name="Dots">
+      <Int Name="Length">96</Int>
+      <Real>0</Real>
+      <Real>0</Real>
+      <Real>-1</Real>
+      <Real>0.18759256558961399</Real>
+      <Real>-0.57735027920788806</Real>
+      <Real>-0.79465444341177649</Real>
+      <Real>-0.49112347245450155</Real>
+      <Real>-0.35682213462759943</Real>
+      <Real>-0.79465445260442358</Real>
+      <Real>0.60706206007773655</Real>
+      <Real>5.2456358414752335e-08</Real>
+      <Real>-0.79465442502648409</Real>
+      <Real>-0.49112352275791088</Real>
+      <Real>0.35682210633554395</Real>
+      <Real>-0.79465443421912996</Real>
+      <Real>0.18759246581169672</Real>
+      <Real>0.57735031162770434</Real>
+      <Real>-0.79465444341177638</Real>
+      <Real>-0.27639310998086691</Real>
+      <Real>-0.85065087452807109</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>0.72360690737127242</Real>
+      <Real>-0.52573101980466264</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>-0.8944272254243314</Real>
+      <Real>0</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>0.72360681651420888</Real>
+      <Real>0.52573114485870043</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>-0.27639325699069922</Real>
+      <Real>0.85065082676168668</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>0.30353113023123213</Real>
+      <Real>-0.93417232306880937</Real>
+      <Real>-0.1875924406599985</Real>
+      <Real>-0.79465443246348733</Real>
+      <Real>-0.57735033372812661</Real>
+      <Real>-0.18759244417127635</Real>
+      <Real>0.98224695410165486</Real>
+      <Real>5.2456354487327486e-08</Real>
+      <Real>-0.18759243363744385</Real>
+      <Real>-0.79465450249337599</Real>
+      <Real>0.57735023962202381</Real>
+      <Real>-0.18759243714872095</Real>
+      <Real>0.30353096878717961</Real>
+      <Real>0.93417237552517096</Real>
+      <Real>-0.18759244065999844</Real>
+      <Real>-0.30353088806515</Real>
+      <Real>-0.93417240175334149</Real>
+      <Real>0.1875924406599985</Real>
+      <Real>0.79465453750831549</Real>
+      <Real>-0.5773501925689698</Real>
+      <Real>0.18759243363744388</Real>
+      <Real>-0.98224695276046647</Real>
+      <Real>-5.245635309796613e-08</Real>
+      <Real>0.18759244065999858</Real>
+      <Real>0.79465446747843327</Real>
+      <Real>0.57735028667507615</Real>
+      <Real>0.18759244065999844</Real>
+      <Real>-0.3035310495092069</Real>
+      <Real>0.93417234929699366</Real>
+      <Real>0.18759244065999844</Real>
+      <Real>0.27639333049561254</Real>
+      <Real>-0.85065080287848482</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>-0.72360677108566906</Real>
+      <Real>-0.52573120738571333</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>0.8944272254243314</Real>
+      <Real>0</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>-0.72360686194274337</Real>
+      <Real>0.52573108233168353</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>0.27639318348578412</Real>
+      <Real>0.85065085064488211</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>-0.18759241592273596</Real>
+      <Real>-0.57735032783760598</Real>
+      <Real>0.79465444341177638</Real>
+      <Real>0.49112354790961382</Real>
+      <Real>-0.35682209218951294</Real>
+      <Real>0.79465442502648398</Real>
+      <Real>-0.60706203601107944</Real>
+      <Real>-5.2456356275319567e-08</Real>
+      <Real>0.79465444341177627</Real>
+      <Real>0.49112349760620688</Real>
+      <Real>0.3568221204815728</Real>
+      <Real>0.79465444341177638</Real>
+      <Real>-0.18759251570065602</Real>
+      <Real>0.57735029541779825</Real>
+      <Real>0.79465444341177627</Real>
+      <Real>0</Real>
+      <Real>0</Real>
+      <Real>1</Real>
+    </Sequence>
+  </SASA>
+</ReferenceData>
diff --git a/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_SurfacePoints42.xml b/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_SurfacePoints42.xml
new file mode 100644 (file)
index 0000000..1f7f0f1
--- /dev/null
@@ -0,0 +1,136 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <SASA Name="Surface">
+    <Real Name="Area">12.566370614359172</Real>
+    <Sequence Name="Dots">
+      <Int Name="Length">126</Int>
+      <Real>0</Real>
+      <Real>0</Real>
+      <Real>-1</Real>
+      <Real>-0.16245979774638578</Real>
+      <Real>-0.5000000507906408</Real>
+      <Real>-0.85065078811787398</Real>
+      <Real>0.42532547872688431</Real>
+      <Real>-0.30901694746435276</Real>
+      <Real>-0.85065078811787398</Real>
+      <Real>-0.52573114485870054</Real>
+      <Real>0</Real>
+      <Real>-0.85065078811787398</Real>
+      <Real>0.42532542532244111</Real>
+      <Real>0.30901702096927364</Real>
+      <Real>-0.85065078811787398</Real>
+      <Real>-0.16245988415659918</Real>
+      <Real>0.5000000227142638</Real>
+      <Real>-0.85065078811787398</Real>
+      <Real>0.26286568098049856</Real>
+      <Real>-0.80901699825499351</Real>
+      <Real>-0.52573104368787149</Real>
+      <Real>-0.68819095823684195</Real>
+      <Real>-0.50000006214777759</Real>
+      <Real>-0.52573105562946909</Real>
+      <Real>0.85065086540550006</Real>
+      <Real>7.3504917541544756e-08</Real>
+      <Real>-0.52573101980467962</Real>
+      <Real>-0.68819101338354338</Real>
+      <Real>0.50000001135712879</Real>
+      <Real>-0.52573103174627511</Real>
+      <Real>0.26286554116584188</Real>
+      <Real>0.80901704368353733</Real>
+      <Real>-0.52573104368787138</Real>
+      <Real>-0.27639310998086691</Real>
+      <Real>-0.85065087452807109</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>0.72360690737127242</Real>
+      <Real>-0.52573101980466264</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>-0.8944272254243314</Real>
+      <Real>0</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>0.72360681651420888</Real>
+      <Real>0.52573114485870043</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>-0.27639325699069922</Real>
+      <Real>0.85065082676168668</Real>
+      <Real>-0.44721352665112024</Real>
+      <Real>1.2961531078992119e-07</Real>
+      <Real>-0.99999999999999167</Real>
+      <Real>5.5511151231257839e-17</Real>
+      <Real>-0.58778516141294568</Real>
+      <Real>-0.80901706040278132</Real>
+      <Real>5.5511151231257827e-17</Real>
+      <Real>0.58778537113492124</Real>
+      <Real>-0.80901690803084136</Real>
+      <Real>-1.1102230246251565e-16</Real>
+      <Real>-0.95105650027381883</Real>
+      <Real>-0.30901704368354138</Real>
+      <Real>1.1102230246251568e-16</Real>
+      <Real>0.95105654032715092</Real>
+      <Real>-0.30901692041205486</Real>
+      <Real>0</Real>
+      <Real>-0.95105652697604204</Real>
+      <Real>0.30901696150255092</Real>
+      <Real>-1.6653345369377348e-16</Real>
+      <Real>0.95105651362493138</Real>
+      <Real>0.30901700259304649</Real>
+      <Real>0</Real>
+      <Real>-0.58778530122760031</Real>
+      <Real>0.80901695882149427</Real>
+      <Real>-1.6653345369377348e-16</Real>
+      <Real>0.58778523132027516</Real>
+      <Real>0.80901700961214085</Real>
+      <Real>0</Real>
+      <Real>-4.3205103550381097e-08</Real>
+      <Real>0.99999999999999911</Real>
+      <Real>0</Real>
+      <Real>0.27639333049561254</Real>
+      <Real>-0.85065080287848482</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>-0.72360677108566906</Real>
+      <Real>-0.52573120738571333</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>0.8944272254243314</Real>
+      <Real>0</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>-0.72360686194274337</Real>
+      <Real>0.52573108233168353</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>0.27639318348578412</Real>
+      <Real>0.85065085064488211</Real>
+      <Real>0.44721352665112024</Real>
+      <Real>-0.26286547125851051</Real>
+      <Real>-0.80901706639780024</Real>
+      <Real>0.52573104368787149</Real>
+      <Real>0.68819104095689176</Real>
+      <Real>-0.49999998596180112</Real>
+      <Real>0.52573101980467962</Real>
+      <Real>-0.85065085064487889</Real>
+      <Real>-7.3504916320299429e-08</Real>
+      <Real>0.52573104368787138</Real>
+      <Real>0.68819098581019356</Real>
+      <Real>0.50000003675245419</Real>
+      <Real>0.52573104368787149</Real>
+      <Real>-0.26286561107317141</Real>
+      <Real>0.80901702096926842</Real>
+      <Real>0.52573104368787149</Real>
+      <Real>0.16245992736170423</Real>
+      <Real>-0.50000000867606953</Real>
+      <Real>0.85065078811787398</Real>
+      <Real>-0.42532539862021468</Real>
+      <Real>-0.30901705772173055</Real>
+      <Real>0.85065078811787398</Real>
+      <Real>0.52573114485870043</Real>
+      <Real>0</Real>
+      <Real>0.85065078811787398</Real>
+      <Real>-0.42532545202466426</Real>
+      <Real>0.30901698421681439</Real>
+      <Real>0.85065078811787398</Real>
+      <Real>0.1624598409514931</Real>
+      <Real>0.50000003675245419</Real>
+      <Real>0.85065078811787398</Real>
+      <Real>0</Real>
+      <Real>0</Real>
+      <Real>1</Real>
+    </Sequence>
+  </SASA>
+</ReferenceData>
index d9eb23d53b062f7d67c041282ca3a3ca3a50920a..8c3fc1745805e7d4745b215c1a1f006e1ebabd7d 100644 (file)
@@ -43,6 +43,8 @@
 
 #include "gromacs/trajectoryanalysis/modules/surfacearea.h"
 
+#include <cstdlib>
+
 #include <gtest/gtest.h>
 
 #include "gromacs/math/utilities.h"
@@ -161,9 +163,20 @@ class SurfaceAreaTest : public ::testing::Test
             dotCount_ = 0;
             sfree(dots_);
             dots_     = NULL;
-            ASSERT_EQ(0, nsc_dclm_pbc(x_, &radius_[0], index_.size(), ndots, flags,
-                                      &area_, &atomArea_, &volume_, &dots_, &dotCount_,
-                                      &index_[0], epbcXYZ, bPBC ? box_ : NULL));
+            t_pbc       pbc;
+            if (bPBC)
+            {
+                set_pbc(&pbc, epbcXYZ, box_);
+            }
+            ASSERT_NO_THROW_GMX(
+                    {
+                        gmx::SurfaceAreaCalculator calculator;
+                        calculator.setDotCount(ndots);
+                        calculator.calculate(x_, &radius_[0], bPBC ? &pbc : NULL,
+                                             index_.size(), &index_[0], flags,
+                                             &area_, &volume_, &atomArea_,
+                                             &dots_, &dotCount_);
+                    });
         }
         real resultArea() const { return area_; }
         real resultVolume() const { return volume_; }
@@ -187,6 +200,10 @@ class SurfaceAreaTest : public ::testing::Test
             {
                 if (checkDotCoordinates)
                 {
+                    // The algorithm may produce the dots in different order in
+                    // single and double precision due to some internal
+                    // sorting...
+                    std::qsort(dots_, dotCount_, sizeof(rvec), &dotComparer);
                     compound.checkSequenceArray(3*dotCount_, dots_, "Dots");
                 }
                 else
@@ -200,6 +217,28 @@ class SurfaceAreaTest : public ::testing::Test
         matrix                          box_;
 
     private:
+        static int dotComparer(const void *a, const void *b)
+        {
+            for (int d = DIM - 1; d >= 0; --d)
+            {
+                const real ad = reinterpret_cast<const real *>(a)[d];
+                const real bd = reinterpret_cast<const real *>(b)[d];
+                // A fudge factor is needed to get an ordering that is the same
+                // in single and double precision, since the points are not
+                // exactly on the same Z plane even though in exact arithmetic
+                // they probably would be.
+                if (ad < bd - 0.001)
+                {
+                    return -1;
+                }
+                else if (ad > bd + 0.001)
+                {
+                    return 1;
+                }
+            }
+            return 0;
+        }
+
         gmx_rng_t          rng_;
         int                allocated_;
         rvec              *x_;
@@ -253,6 +292,42 @@ TEST_F(SurfaceAreaTest, ComputesTwoPointsOfUnequalRadius)
     EXPECT_REAL_EQ_TOL(2*M_PI*(dist - 0.5 + 2)*2, atomArea(1), tolerance);
 }
 
+TEST_F(SurfaceAreaTest, SurfacePoints12)
+{
+    gmx::test::TestReferenceChecker checker(data_.rootChecker());
+    reserveSpace(1);
+    addSphere(0, 0, 0, 1);
+    ASSERT_NO_FATAL_FAILURE(calculate(12, FLAG_DOTS, false));
+    checkReference(&checker, "Surface", true);
+}
+
+TEST_F(SurfaceAreaTest, SurfacePoints32)
+{
+    gmx::test::TestReferenceChecker checker(data_.rootChecker());
+    reserveSpace(1);
+    addSphere(0, 0, 0, 1);
+    ASSERT_NO_FATAL_FAILURE(calculate(32, FLAG_DOTS, false));
+    checkReference(&checker, "Surface", true);
+}
+
+TEST_F(SurfaceAreaTest, SurfacePoints42)
+{
+    gmx::test::TestReferenceChecker checker(data_.rootChecker());
+    reserveSpace(1);
+    addSphere(0, 0, 0, 1);
+    ASSERT_NO_FATAL_FAILURE(calculate(42, FLAG_DOTS, false));
+    checkReference(&checker, "Surface", true);
+}
+
+TEST_F(SurfaceAreaTest, SurfacePoints122)
+{
+    gmx::test::TestReferenceChecker checker(data_.rootChecker());
+    reserveSpace(1);
+    addSphere(0, 0, 0, 1);
+    ASSERT_NO_FATAL_FAILURE(calculate(122, FLAG_DOTS, false));
+    checkReference(&checker, "Surface", true);
+}
+
 TEST_F(SurfaceAreaTest, Computes100Points)
 {
     gmx::test::TestReferenceChecker checker(data_.rootChecker());