Add some unit tests for surface area calculation

There are a few real unit tests that compare the results against
analytical calculation for trivial cases, but the rest of the tests are
more of regression tests, designed to catch changes in the output.
I have not actually checked that the reference values for the more
complex values are correct.

Put the tests in surfacearea.cpp instead of nsc.cpp and renamed the
implementation to match in preparation for future cleanup.

The volume computation with PBC is not tested, since it is broken.

Change-Id: I18e3baa899b1c10fec0c8b43d8dd76a357297005

diff --git a/src/gromacs/trajectoryanalysis/modules/sasa.cpp b/src/gromacs/trajectoryanalysis/modules/sasa.cpp
index 3e5f2543d813187a171f2b91daaa666dbfd03765..70654e76c9b1d306bd8d7bffd03fb47f61a73a8f 100644 (file)
--- a/src/gromacs/trajectoryanalysis/modules/sasa.cpp
+++ b/src/gromacs/trajectoryanalysis/modules/sasa.cpp
@@ -75,7 +75,7 @@
 #include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/stringutil.h"
 
-#include "nsc.h"
+#include "surfacearea.h"
 
 namespace gmx
 {

diff --git a/src/gromacs/trajectoryanalysis/modules/nsc.c b/src/gromacs/trajectoryanalysis/modules/surfacearea.c
similarity index 99%
rename from src/gromacs/trajectoryanalysis/modules/nsc.c
rename to src/gromacs/trajectoryanalysis/modules/surfacearea.c
index 29f9875b6e10e28c4d213dfb7f176ea40c13e0ff..1c6c6f2a73092d32f21e8f425a1534d591c988d8 100644 (file)
--- a/src/gromacs/trajectoryanalysis/modules/nsc.c
+++ b/src/gromacs/trajectoryanalysis/modules/surfacearea.c
@@ -36,7 +36,7 @@
  */
 #include "gmxpre.h"
 
-#include "nsc.h"
+#include "surfacearea.h"
 
 #include <math.h>
 #include <stdarg.h>

diff --git a/src/gromacs/trajectoryanalysis/modules/nsc.h b/src/gromacs/trajectoryanalysis/modules/surfacearea.h
similarity index 96%
rename from src/gromacs/trajectoryanalysis/modules/nsc.h
rename to src/gromacs/trajectoryanalysis/modules/surfacearea.h
index bb39b9f35053f8a8ee2eefa9a8fd41553f091905..eb5129ea22d8de6d0da3ea67c35a676a78a02371 100644 (file)
--- a/src/gromacs/trajectoryanalysis/modules/nsc.h
+++ b/src/gromacs/trajectoryanalysis/modules/surfacearea.h
@@ -34,8 +34,8 @@
  * To help us fund GROMACS development, we humbly ask that you cite
  * the research papers on the package. Check out http://www.gromacs.org.
  */
-#ifndef GMX_TRAJECTORYANALYSIS_NSC_H
-#define GMX_TRAJECTORYANALYSIS_NSC_H
+#ifndef GMX_TRAJECTORYANALYSIS_SURFACEAREA_H
+#define GMX_TRAJECTORYANALYSIS_SURFACEAREA_H
 
 #include "gromacs/legacyheaders/types/simple.h"
 
@@ -107,13 +107,13 @@ int nsc_dclm_pbc(const rvec *coords, real *radius, int nat,
    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 nsc.h !!
+   All files calling NSC should include surfacearea.h !!
 
 
    Example for calling NSC (contents of user file):
 
    ...
-   #include "nsc.h"
+   #include "surfacearea.h"
 
    int routine_calling_NSC(int n_atom, real *coordinates, real *radii) {
    real area, volume, *atomwise_area, *surface_dots;

diff --git a/src/gromacs/trajectoryanalysis/tests/CMakeLists.txt b/src/gromacs/trajectoryanalysis/tests/CMakeLists.txt
index b07003f628bef568c0266ca6f94a370261fd52a4..9d00ed3300ee7c9249aeb760b74a2e01ffc9531b 100644 (file)
--- a/src/gromacs/trajectoryanalysis/tests/CMakeLists.txt
+++ b/src/gromacs/trajectoryanalysis/tests/CMakeLists.txt
@@ -39,6 +39,7 @@ gmx_add_unit_test(TrajectoryAnalysisUnitTests trajectoryanalysis-test
                   freevolume.cpp
                   sasa.cpp
                   select.cpp
+                  surfacearea.cpp
                   $<TARGET_OBJECTS:analysisdata-test-shared>)
 
 add_executable(test_selection ${UNITTEST_TARGET_OPTIONS} test_selection.cpp)

diff --git a/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_Computes100Points.xml b/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_Computes100Points.xml
new file mode 100644 (file)
index 0000000..490019c
--- /dev/null
+++ b/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_Computes100Points.xml
@@ -0,0 +1,112 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <SASA Name="100Points">
+    <Real Name="Area">1041.9522794772724</Real>
+    <Real Name="Volume">808.03995831604618</Real>
+    <Sequence Name="AtomArea">
+      <Int Name="Length">100</Int>
+      <Real>0.95131516208197198</Real>
+      <Real>19.01541145000165</Real>
+      <Real>22.043110642171214</Real>
+      <Real>18.151215320405854</Real>
+      <Real>0.63204646296068279</Real>
+      <Real>0</Real>
+      <Real>15.153269443643612</Real>
+      <Real>7.3352383071855414</Real>
+      <Real>4.1967272329122345</Real>
+      <Real>0</Real>
+      <Real>14.120138371629922</Real>
+      <Real>17.613045388151775</Real>
+      <Real>24.026150792470702</Real>
+      <Real>18.227081273770153</Real>
+      <Real>3.4343773695403952</Real>
+      <Real>19.754449542129247</Real>
+      <Real>21.485580504254141</Real>
+      <Real>12.292225741625185</Real>
+      <Real>19.030725859555776</Real>
+      <Real>0</Real>
+      <Real>9.4242992760866073</Real>
+      <Real>2.3664904848604777</Real>
+      <Real>21.704207650824099</Real>
+      <Real>0</Real>
+      <Real>13.730659331066482</Real>
+      <Real>14.531937785202674</Real>
+      <Real>4.5894783541023525</Real>
+      <Real>3.6755820832244703</Real>
+      <Real>0.44427283053350852</Real>
+      <Real>0.1066189700313567</Real>
+      <Real>2.0456369215234842</Real>
+      <Real>8.5740365276665162</Real>
+      <Real>18.165852189840816</Real>
+      <Real>10.496616250155016</Real>
+      <Real>21.080864118880193</Real>
+      <Real>2.3899727797298986</Real>
+      <Real>0.78849936521174724</Real>
+      <Real>13.408721074833956</Real>
+      <Real>13.658708645471602</Real>
+      <Real>25.966986031945716</Real>
+      <Real>0.85119617060510389</Real>
+      <Real>19.63570848962739</Real>
+      <Real>0.14968475148186061</Real>
+      <Real>10.306327915331671</Real>
+      <Real>4.3963272956306936</Real>
+      <Real>24.525460566049077</Real>
+      <Real>8.0156199763235296</Real>
+      <Real>0</Real>
+      <Real>5.3750282656324613</Real>
+      <Real>0</Real>
+      <Real>12.294379373044716</Real>
+      <Real>10.315417503550229</Real>
+      <Real>12.029185950445784</Real>
+      <Real>2.0160888710349627</Real>
+      <Real>0</Real>
+      <Real>0.13787251489085933</Real>
+      <Real>7.9417939895007734</Real>
+      <Real>42.12368047628528</Real>
+      <Real>15.724000321454668</Real>
+      <Real>12.462533452088378</Real>
+      <Real>15.79912862220181</Real>
+      <Real>6.2688853022092417</Real>
+      <Real>12.486923176542316</Real>
+      <Real>6.6320474695863378</Real>
+      <Real>20.983750899642256</Real>
+      <Real>30.213450121453786</Real>
+      <Real>0.6401170350458556</Real>
+      <Real>5.4877256580120797</Real>
+      <Real>7.5123761006137872</Real>
+      <Real>6.6756067906869934</Real>
+      <Real>14.293717932888253</Real>
+      <Real>6.9724919386790729</Real>
+      <Real>27.362816599110033</Real>
+      <Real>4.4291554575110448</Real>
+      <Real>3.7993133931974961</Real>
+      <Real>8.2368423448648969</Real>
+      <Real>6.5443862216435997</Real>
+      <Real>5.4944959373624975</Real>
+      <Real>31.691642853696951</Real>
+      <Real>12.180941103950339</Real>
+      <Real>24.133979864235052</Real>
+      <Real>11.141753032303363</Real>
+      <Real>12.053564446854951</Real>
+      <Real>3.0590924364990446</Real>
+      <Real>5.0901076005888326</Real>
+      <Real>12.701834505720841</Real>
+      <Real>1.6067117447360322</Real>
+      <Real>8.586295906028095</Real>
+      <Real>3.7715811132613966</Real>
+      <Real>13.891761175025035</Real>
+      <Real>14.091748610057952</Real>
+      <Real>13.276603955659391</Real>
+      <Real>9.7895079276620454</Real>
+      <Real>0.80200851184333233</Real>
+      <Real>0</Real>
+      <Real>17.970129415135109</Real>
+      <Real>24.853015014087575</Real>
+      <Real>4.7620404205235838</Real>
+      <Real>9.4941606968585042</Real>
+      <Real>4.2587127226353019</Real>
+    </Sequence>
+    <Int Name="DotCount">1338</Int>
+  </SASA>
+</ReferenceData>

diff --git a/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_Computes100PointsWithRectangularPBC.xml b/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_Computes100PointsWithRectangularPBC.xml
new file mode 100644 (file)
index 0000000..c30bb64
--- /dev/null
+++ b/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_Computes100PointsWithRectangularPBC.xml
@@ -0,0 +1,111 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <SASA Name="100Points">
+    <Real Name="Area">1041.9522794772727</Real>
+    <Sequence Name="AtomArea">
+      <Int Name="Length">100</Int>
+      <Real>0.95131516208197198</Real>
+      <Real>19.01541145000165</Real>
+      <Real>22.043110642171214</Real>
+      <Real>18.151215320405854</Real>
+      <Real>0.63204646296068279</Real>
+      <Real>0</Real>
+      <Real>15.153269443643612</Real>
+      <Real>7.3352383071855414</Real>
+      <Real>4.1967272329122345</Real>
+      <Real>0</Real>
+      <Real>14.120138371629922</Real>
+      <Real>17.613045388151775</Real>
+      <Real>24.026150792470702</Real>
+      <Real>18.227081273770153</Real>
+      <Real>3.4343773695403952</Real>
+      <Real>19.754449542129247</Real>
+      <Real>21.485580504254141</Real>
+      <Real>12.292225741625185</Real>
+      <Real>19.030725859555776</Real>
+      <Real>0</Real>
+      <Real>9.4242992760866073</Real>
+      <Real>2.3664904848604777</Real>
+      <Real>21.704207650824099</Real>
+      <Real>0</Real>
+      <Real>13.730659331066482</Real>
+      <Real>14.531937785202674</Real>
+      <Real>4.5894783541023525</Real>
+      <Real>3.6755820832244703</Real>
+      <Real>0.44427283053350852</Real>
+      <Real>0.1066189700313567</Real>
+      <Real>2.0456369215234842</Real>
+      <Real>8.5740365276665162</Real>
+      <Real>18.165852189840816</Real>
+      <Real>10.496616250155016</Real>
+      <Real>21.080864118880193</Real>
+      <Real>2.3899727797298986</Real>
+      <Real>0.78849936521174724</Real>
+      <Real>13.408721074833956</Real>
+      <Real>13.658708645471602</Real>
+      <Real>25.966986031945716</Real>
+      <Real>0.85119617060510389</Real>
+      <Real>19.63570848962739</Real>
+      <Real>0.14968475148186061</Real>
+      <Real>10.306327915331671</Real>
+      <Real>4.3963272956306936</Real>
+      <Real>24.525460566049077</Real>
+      <Real>8.0156199763235296</Real>
+      <Real>0</Real>
+      <Real>5.3750282656324613</Real>
+      <Real>0</Real>
+      <Real>12.294379373044716</Real>
+      <Real>10.315417503550229</Real>
+      <Real>12.029185950445784</Real>
+      <Real>2.0160888710349627</Real>
+      <Real>0</Real>
+      <Real>0.13787251489085933</Real>
+      <Real>7.9417939895007734</Real>
+      <Real>42.12368047628528</Real>
+      <Real>15.724000321454668</Real>
+      <Real>12.462533452088378</Real>
+      <Real>15.79912862220181</Real>
+      <Real>6.2688853022092417</Real>
+      <Real>12.486923176542316</Real>
+      <Real>6.6320474695863378</Real>
+      <Real>20.983750899642256</Real>
+      <Real>30.213450121453786</Real>
+      <Real>0.6401170350458556</Real>
+      <Real>5.4877256580120797</Real>
+      <Real>7.5123761006137872</Real>
+      <Real>6.6756067906869934</Real>
+      <Real>14.293717932888253</Real>
+      <Real>6.9724919386790729</Real>
+      <Real>27.362816599110033</Real>
+      <Real>4.4291554575110448</Real>
+      <Real>3.7993133931974961</Real>
+      <Real>8.2368423448648969</Real>
+      <Real>6.5443862216435997</Real>
+      <Real>5.4944959373624975</Real>
+      <Real>31.691642853696951</Real>
+      <Real>12.180941103950339</Real>
+      <Real>24.133979864235052</Real>
+      <Real>11.141753032303363</Real>
+      <Real>12.053564446854951</Real>
+      <Real>3.0590924364990446</Real>
+      <Real>5.0901076005888326</Real>
+      <Real>12.701834505720841</Real>
+      <Real>1.6067117447360322</Real>
+      <Real>8.586295906028095</Real>
+      <Real>3.7715811132613966</Real>
+      <Real>13.891761175025035</Real>
+      <Real>14.091748610057952</Real>
+      <Real>13.276603955659391</Real>
+      <Real>9.7895079276620454</Real>
+      <Real>0.80200851184333233</Real>
+      <Real>0</Real>
+      <Real>17.970129415135109</Real>
+      <Real>24.853015014087575</Real>
+      <Real>4.7620404205235838</Real>
+      <Real>9.4941606968585042</Real>
+      <Real>4.2587127226353019</Real>
+    </Sequence>
+    <Int Name="DotCount">1338</Int>
+  </SASA>
+</ReferenceData>

diff --git a/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_Computes100PointsWithTriclinicPBC.xml b/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_Computes100PointsWithTriclinicPBC.xml
new file mode 100644 (file)
index 0000000..c30bb64
--- /dev/null
+++ b/src/gromacs/trajectoryanalysis/tests/refdata/SurfaceAreaTest_Computes100PointsWithTriclinicPBC.xml
@@ -0,0 +1,111 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <SASA Name="100Points">
+    <Real Name="Area">1041.9522794772727</Real>
+    <Sequence Name="AtomArea">
+      <Int Name="Length">100</Int>
+      <Real>0.95131516208197198</Real>
+      <Real>19.01541145000165</Real>
+      <Real>22.043110642171214</Real>
+      <Real>18.151215320405854</Real>
+      <Real>0.63204646296068279</Real>
+      <Real>0</Real>
+      <Real>15.153269443643612</Real>
+      <Real>7.3352383071855414</Real>
+      <Real>4.1967272329122345</Real>
+      <Real>0</Real>
+      <Real>14.120138371629922</Real>
+      <Real>17.613045388151775</Real>
+      <Real>24.026150792470702</Real>
+      <Real>18.227081273770153</Real>
+      <Real>3.4343773695403952</Real>
+      <Real>19.754449542129247</Real>
+      <Real>21.485580504254141</Real>
+      <Real>12.292225741625185</Real>
+      <Real>19.030725859555776</Real>
+      <Real>0</Real>
+      <Real>9.4242992760866073</Real>
+      <Real>2.3664904848604777</Real>
+      <Real>21.704207650824099</Real>
+      <Real>0</Real>
+      <Real>13.730659331066482</Real>
+      <Real>14.531937785202674</Real>
+      <Real>4.5894783541023525</Real>
+      <Real>3.6755820832244703</Real>
+      <Real>0.44427283053350852</Real>
+      <Real>0.1066189700313567</Real>
+      <Real>2.0456369215234842</Real>
+      <Real>8.5740365276665162</Real>
+      <Real>18.165852189840816</Real>
+      <Real>10.496616250155016</Real>
+      <Real>21.080864118880193</Real>
+      <Real>2.3899727797298986</Real>
+      <Real>0.78849936521174724</Real>
+      <Real>13.408721074833956</Real>
+      <Real>13.658708645471602</Real>
+      <Real>25.966986031945716</Real>
+      <Real>0.85119617060510389</Real>
+      <Real>19.63570848962739</Real>
+      <Real>0.14968475148186061</Real>
+      <Real>10.306327915331671</Real>
+      <Real>4.3963272956306936</Real>
+      <Real>24.525460566049077</Real>
+      <Real>8.0156199763235296</Real>
+      <Real>0</Real>
+      <Real>5.3750282656324613</Real>
+      <Real>0</Real>
+      <Real>12.294379373044716</Real>
+      <Real>10.315417503550229</Real>
+      <Real>12.029185950445784</Real>
+      <Real>2.0160888710349627</Real>
+      <Real>0</Real>
+      <Real>0.13787251489085933</Real>
+      <Real>7.9417939895007734</Real>
+      <Real>42.12368047628528</Real>
+      <Real>15.724000321454668</Real>
+      <Real>12.462533452088378</Real>
+      <Real>15.79912862220181</Real>
+      <Real>6.2688853022092417</Real>
+      <Real>12.486923176542316</Real>
+      <Real>6.6320474695863378</Real>
+      <Real>20.983750899642256</Real>
+      <Real>30.213450121453786</Real>
+      <Real>0.6401170350458556</Real>
+      <Real>5.4877256580120797</Real>
+      <Real>7.5123761006137872</Real>
+      <Real>6.6756067906869934</Real>
+      <Real>14.293717932888253</Real>
+      <Real>6.9724919386790729</Real>
+      <Real>27.362816599110033</Real>
+      <Real>4.4291554575110448</Real>
+      <Real>3.7993133931974961</Real>
+      <Real>8.2368423448648969</Real>
+      <Real>6.5443862216435997</Real>
+      <Real>5.4944959373624975</Real>
+      <Real>31.691642853696951</Real>
+      <Real>12.180941103950339</Real>
+      <Real>24.133979864235052</Real>
+      <Real>11.141753032303363</Real>
+      <Real>12.053564446854951</Real>
+      <Real>3.0590924364990446</Real>
+      <Real>5.0901076005888326</Real>
+      <Real>12.701834505720841</Real>
+      <Real>1.6067117447360322</Real>
+      <Real>8.586295906028095</Real>
+      <Real>3.7715811132613966</Real>
+      <Real>13.891761175025035</Real>
+      <Real>14.091748610057952</Real>
+      <Real>13.276603955659391</Real>
+      <Real>9.7895079276620454</Real>
+      <Real>0.80200851184333233</Real>
+      <Real>0</Real>
+      <Real>17.970129415135109</Real>
+      <Real>24.853015014087575</Real>
+      <Real>4.7620404205235838</Real>
+      <Real>9.4941606968585042</Real>
+      <Real>4.2587127226353019</Real>
+    </Sequence>
+    <Int Name="DotCount">1338</Int>
+  </SASA>
+</ReferenceData>

diff --git a/src/gromacs/trajectoryanalysis/tests/sasa.cpp b/src/gromacs/trajectoryanalysis/tests/sasa.cpp
index 75a60d69df988d07239464b25129885138fec9f5..55515b411568350738622349ff21b4da0d0323d4 100644 (file)
--- a/src/gromacs/trajectoryanalysis/tests/sasa.cpp
+++ b/src/gromacs/trajectoryanalysis/tests/sasa.cpp
@@ -46,8 +46,7 @@
  *  - Tests for XVG labels.  This is a limitation of the current testing
  *    framework.
  *
- * The actual surface area algorithm should be tested separately with a more
- * extensive set of data, but those fit better as a separate set of unit tests.
+ * The actual surface area algorithm is tested separately in surfacearea.cpp.
  *
  * \author Teemu Murtola <teemu.murtola@gmail.com>
  * \ingroup module_trajectoryanalysis

diff --git a/src/gromacs/trajectoryanalysis/tests/surfacearea.cpp b/src/gromacs/trajectoryanalysis/tests/surfacearea.cpp
new file mode 100644 (file)
index 0000000..d9eb23d
--- /dev/null
+++ b/src/gromacs/trajectoryanalysis/tests/surfacearea.cpp
@@ -0,0 +1,334 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Tests for the surface area calculation used by the `sasa` analysis module.
+ *
+ * \author Teemu Murtola <teemu.murtola@gmail.com>
+ * \ingroup module_trajectoryanalysis
+ */
+#include "gmxpre.h"
+
+#include "gromacs/trajectoryanalysis/modules/surfacearea.h"
+
+#include <gtest/gtest.h>
+
+#include "gromacs/math/utilities.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/random/random.h"
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/smalloc.h"
+
+#include "testutils/refdata.h"
+#include "testutils/testasserts.h"
+
+namespace
+{
+
+/********************************************************************
+ * SurfaceAreaTest
+ */
+
+class SurfaceAreaTest : public ::testing::Test
+{
+    public:
+        SurfaceAreaTest()
+            : rng_(NULL), allocated_(0), x_(NULL), area_(0.0), volume_(0.0),
+              atomArea_(NULL), dotCount_(0), dots_(NULL)
+        {
+            // TODO: Handle errors.
+            rng_ = gmx_rng_init(12345);
+            clear_mat(box_);
+        }
+        ~SurfaceAreaTest()
+        {
+            if (rng_ != NULL)
+            {
+                gmx_rng_destroy(rng_);
+            }
+            sfree(x_);
+            sfree(atomArea_);
+            sfree(dots_);
+        }
+
+        void reserveSpace(int count)
+        {
+            GMX_RELEASE_ASSERT(allocated_ == 0 && x_ == NULL,
+                               "Cannot allocate data more than once");
+            snew(x_, count);
+            radius_.reserve(count);
+            index_.reserve(count);
+            allocated_ = count;
+        }
+
+        void addSphere(real x, real y, real z, real radius,
+                       bool bAddToIndex = true)
+        {
+            const int index = radius_.size();
+            GMX_RELEASE_ASSERT(index < allocated_,
+                               "Not enough space allocated for test data");
+            x_[index][XX] = x;
+            x_[index][YY] = y;
+            x_[index][ZZ] = z;
+            if (bAddToIndex)
+            {
+                index_.push_back(radius_.size());
+            }
+            radius_.push_back(radius);
+        }
+
+        void generateRandomPosition(rvec x, real *radius)
+        {
+            rvec fx;
+            fx[XX]  = gmx_rng_uniform_real(rng_);
+            fx[YY]  = gmx_rng_uniform_real(rng_);
+            fx[ZZ]  = gmx_rng_uniform_real(rng_);
+            mvmul(box_, fx, x);
+            *radius = 1.5*gmx_rng_uniform_real(rng_) + 0.5;
+        }
+
+        void addDummySpheres(int count)
+        {
+            for (int i = 0; i < count; ++i)
+            {
+                rvec x;
+                real radius;
+                generateRandomPosition(x, &radius);
+                addSphere(x[XX], x[YY], x[ZZ], radius, false);
+            }
+        }
+
+        void generateRandomPositions(int count)
+        {
+            reserveSpace(count);
+            for (int i = 0; i < count; ++i)
+            {
+                rvec x;
+                real radius;
+                generateRandomPosition(x, &radius);
+                addSphere(x[XX], x[YY], x[ZZ], radius);
+            }
+        }
+        void translatePoints(real x, real y, real z)
+        {
+            for (size_t i = 0; i < radius_.size(); ++i)
+            {
+                x_[i][XX] += x;
+                x_[i][YY] += y;
+                x_[i][ZZ] += z;
+            }
+        }
+
+        void calculate(int ndots, int flags, bool bPBC)
+        {
+            volume_   = 0.0;
+            sfree(atomArea_);
+            atomArea_ = NULL;
+            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));
+        }
+        real resultArea() const { return area_; }
+        real resultVolume() const { return volume_; }
+        real atomArea(int index) const { return atomArea_[index]; }
+
+        void checkReference(gmx::test::TestReferenceChecker *checker, const char *id,
+                            bool checkDotCoordinates)
+        {
+            gmx::test::TestReferenceChecker compound(
+                    checker->checkCompound("SASA", id));
+            compound.checkReal(area_, "Area");
+            if (volume_ > 0.0)
+            {
+                compound.checkReal(volume_, "Volume");
+            }
+            if (atomArea_ != NULL)
+            {
+                compound.checkSequenceArray(index_.size(), atomArea_, "AtomArea");
+            }
+            if (dots_ != NULL)
+            {
+                if (checkDotCoordinates)
+                {
+                    compound.checkSequenceArray(3*dotCount_, dots_, "Dots");
+                }
+                else
+                {
+                    compound.checkInteger(dotCount_, "DotCount");
+                }
+            }
+        }
+
+        gmx::test::TestReferenceData    data_;
+        matrix                          box_;
+
+    private:
+        gmx_rng_t          rng_;
+        int                allocated_;
+        rvec              *x_;
+        std::vector<real>  radius_;
+        std::vector<int>   index_;
+
+        real               area_;
+        real               volume_;
+        real              *atomArea_;
+        int                dotCount_;
+        real              *dots_;
+};
+
+TEST_F(SurfaceAreaTest, ComputesSinglePoint)
+{
+    gmx::test::FloatingPointTolerance tolerance(
+            gmx::test::defaultRealTolerance());
+    reserveSpace(1);
+    addSphere(1, 1, 1, 1);
+    ASSERT_NO_FATAL_FAILURE(calculate(24, FLAG_VOLUME | FLAG_ATOM_AREA, false));
+    EXPECT_REAL_EQ_TOL(4*M_PI, resultArea(), tolerance);
+    EXPECT_REAL_EQ_TOL(4*M_PI, atomArea(0), tolerance);
+    EXPECT_REAL_EQ_TOL(4*M_PI/3, resultVolume(), tolerance);
+}
+
+TEST_F(SurfaceAreaTest, ComputesTwoPoints)
+{
+    gmx::test::FloatingPointTolerance tolerance(
+            gmx::test::relativeToleranceAsFloatingPoint(1.0, 0.005));
+    reserveSpace(2);
+    addSphere(1, 1, 1, 1);
+    addSphere(2, 1, 1, 1);
+    ASSERT_NO_FATAL_FAILURE(calculate(1000, FLAG_ATOM_AREA, false));
+    EXPECT_REAL_EQ_TOL(2*2*M_PI*1.5, resultArea(), tolerance);
+    EXPECT_REAL_EQ_TOL(2*M_PI*1.5, atomArea(0), tolerance);
+    EXPECT_REAL_EQ_TOL(2*M_PI*1.5, atomArea(1), tolerance);
+}
+
+TEST_F(SurfaceAreaTest, ComputesTwoPointsOfUnequalRadius)
+{
+    gmx::test::FloatingPointTolerance tolerance(
+            gmx::test::relativeToleranceAsFloatingPoint(1.0, 0.005));
+    reserveSpace(2);
+    // Spheres of radius 1 and 2 with intersection at 1.5
+    const real dist = 0.5 + sqrt(3.25);
+    addSphere(1.0, 1.0, 1.0, 1);
+    addSphere(1.0 + dist, 1.0, 1.0, 2);
+    ASSERT_NO_FATAL_FAILURE(calculate(1000, FLAG_ATOM_AREA, false));
+    EXPECT_REAL_EQ_TOL(2*M_PI*(1.5 + (dist - 0.5 + 2)*2), resultArea(), tolerance);
+    EXPECT_REAL_EQ_TOL(2*M_PI*1.5, atomArea(0), tolerance);
+    EXPECT_REAL_EQ_TOL(2*M_PI*(dist - 0.5 + 2)*2, atomArea(1), tolerance);
+}
+
+TEST_F(SurfaceAreaTest, Computes100Points)
+{
+    gmx::test::TestReferenceChecker checker(data_.rootChecker());
+    checker.setDefaultTolerance(gmx::test::absoluteTolerance(0.001));
+    box_[XX][XX] = 10.0;
+    box_[YY][YY] = 10.0;
+    box_[ZZ][ZZ] = 10.0;
+    generateRandomPositions(100);
+    ASSERT_NO_FATAL_FAILURE(calculate(24, FLAG_VOLUME | FLAG_ATOM_AREA | FLAG_DOTS, false));
+    checkReference(&checker, "100Points", false);
+}
+
+TEST_F(SurfaceAreaTest, Computes100PointsWithRectangularPBC)
+{
+    // TODO: It would be nice to check that this produces the same result as
+    // without PBC, without duplicating the reference files.
+    gmx::test::TestReferenceChecker checker(data_.rootChecker());
+    checker.setDefaultTolerance(gmx::test::absoluteTolerance(0.001));
+    box_[XX][XX] = 10.0;
+    box_[YY][YY] = 10.0;
+    box_[ZZ][ZZ] = 10.0;
+    generateRandomPositions(100);
+    box_[XX][XX] = 20.0;
+    box_[YY][YY] = 20.0;
+    box_[ZZ][ZZ] = 20.0;
+    // TODO: The volume computation with PBC is broken...
+    const int flags = FLAG_ATOM_AREA | FLAG_DOTS;
+    ASSERT_NO_FATAL_FAILURE(calculate(24, flags, true));
+    checkReference(&checker, "100Points", false);
+
+    translatePoints(15.0, 0, 0);
+    ASSERT_NO_FATAL_FAILURE(calculate(24, flags, true));
+    checkReference(&checker, "100Points", false);
+
+    translatePoints(-15.0, 15.0, 0);
+    ASSERT_NO_FATAL_FAILURE(calculate(24, flags, true));
+    checkReference(&checker, "100Points", false);
+
+    translatePoints(0, -15.0, 15.0);
+    ASSERT_NO_FATAL_FAILURE(calculate(24, flags, true));
+    checkReference(&checker, "100Points", false);
+}
+
+TEST_F(SurfaceAreaTest, Computes100PointsWithTriclinicPBC)
+{
+    // TODO: It would be nice to check that this produces the same result as
+    // without PBC, without duplicating the reference files.
+    gmx::test::TestReferenceChecker checker(data_.rootChecker());
+    checker.setDefaultTolerance(gmx::test::absoluteTolerance(0.001));
+    box_[XX][XX] = 10.0;
+    box_[YY][YY] = 10.0;
+    box_[ZZ][ZZ] = 10.0;
+    generateRandomPositions(100);
+    box_[XX][XX] = 20.0;
+    box_[YY][XX] = 10.0;
+    box_[YY][YY] = 10.0*sqrt(3.0);
+    box_[ZZ][XX] = 10.0;
+    box_[ZZ][YY] = 10.0*sqrt(1.0/3.0);
+    box_[ZZ][ZZ] = 20.0*sqrt(2.0/3.0);
+
+    // TODO: The volume computation with PBC is broken...
+    const int flags = FLAG_ATOM_AREA | FLAG_DOTS;
+    ASSERT_NO_FATAL_FAILURE(calculate(24, flags, true));
+    checkReference(&checker, "100Points", false);
+
+    translatePoints(15.0, 0, 0);
+    ASSERT_NO_FATAL_FAILURE(calculate(24, flags, true));
+    checkReference(&checker, "100Points", false);
+
+    translatePoints(-15.0, box_[YY][YY] - 5.0, 0);
+    ASSERT_NO_FATAL_FAILURE(calculate(24, flags, true));
+    checkReference(&checker, "100Points", false);
+
+    translatePoints(0, -(box_[YY][YY] - 5.0), 15.0);
+    ASSERT_NO_FATAL_FAILURE(calculate(24, flags, true));
+    checkReference(&checker, "100Points", false);
+}
+
+} // namespace

diff --git a/src/testutils/testasserts.h b/src/testutils/testasserts.h
index e38600fb8b0ba02e763d80edd7dec9e55de80882..bce7395e5249a4e7655e6cda4bc0b107870d94bc 100644 (file)
--- a/src/testutils/testasserts.h
+++ b/src/testutils/testasserts.h
@@ -401,6 +401,18 @@ relativeToleranceAsPrecisionDependentUlp(double       magnitude,
                                   singleUlpDiff, doubleUlpDiff, false);
 }
 
+/*! \brief
+ * Creates a tolerance that allows a specified absolute difference.
+ *
+ * \related FloatingPointTolerance
+ */
+static inline FloatingPointTolerance
+absoluteTolerance(double tolerance)
+{
+    return FloatingPointTolerance(tolerance, tolerance,
+                                  GMX_UINT64_MAX, GMX_UINT64_MAX, false);
+}
+
 /*! \brief
  * Creates a tolerance that allows a relative difference in a complex
  * computation.