Convert gmx rdf to the C++ framework
authorTeemu Murtola <teemu.murtola@gmail.com>
Thu, 26 Jun 2014 11:59:07 +0000 (14:59 +0300)
committerGerrit Code Review <gerrit@gerrit.gromacs.org>
Thu, 11 Dec 2014 12:29:19 +0000 (13:29 +0100)
Convert the RDF calculation tool to the new framework.  Main changed
beyond mechanical conversion:
 - Remove -com and -rdf options in favor of using selections directly.
   Both functionality, and a lot more, is now directly supported with
   proper selections.
 - Remove -hq and -fade options, as both are just post-processing that
   can be trivially done based on the output xvg file alone.
 - Normalization is reworked to work with dynamic selections.
 - Add an -rmax parameter that allows specifying an upper cutoff for the
   RDF.  Combined with the use of the neighborhood search routines, this
   potentially makes the tool a lot faster if the RDF is not required up
   to half the box length.

Add basic unit/regression tests for a simple water box.  The reference
values are not verified in any way, so they mainly catch regression
bugs.

Change-Id: Ieed19fc39a84e299291232190333652ced26abee

16 files changed:
src/gromacs/analysisdata/modules/plot.cpp
src/gromacs/analysisdata/modules/plot.h
src/gromacs/gmxana/gmx_ana.h
src/gromacs/gmxana/gmx_rdf.c [deleted file]
src/gromacs/trajectoryanalysis/modules.cpp
src/gromacs/trajectoryanalysis/modules/rdf.cpp [new file with mode: 0644]
src/gromacs/trajectoryanalysis/modules/rdf.h [new file with mode: 0644]
src/gromacs/trajectoryanalysis/tests/CMakeLists.txt
src/gromacs/trajectoryanalysis/tests/rdf.cpp [new file with mode: 0644]
src/gromacs/trajectoryanalysis/tests/refdata/RdfModuleTest_BasicTest.xml [new file with mode: 0644]
src/gromacs/trajectoryanalysis/tests/refdata/RdfModuleTest_CalculatesSurf.xml [new file with mode: 0644]
src/gromacs/trajectoryanalysis/tests/refdata/RdfModuleTest_CalculatesXY.xml [new file with mode: 0644]
src/gromacs/trajectoryanalysis/tests/sasa.cpp
src/gromacs/trajectoryanalysis/tests/spc216.gro [new file with mode: 0644]
src/programs/CreateLinks.cmake.cmakein
src/programs/legacymodules.cpp

index e421e2d64a85aea39d64679b758993531f6265a9..e60648e4fa333b6275fca1396b92755817598e3a 100644 (file)
@@ -226,6 +226,12 @@ AbstractPlotModule::setTitle(const char *title)
     impl_->title_ = title;
 }
 
+void
+AbstractPlotModule::setTitle(const std::string &title)
+{
+    impl_->title_ = title;
+}
+
 
 void
 AbstractPlotModule::setSubtitle(const char *subtitle)
@@ -234,6 +240,13 @@ AbstractPlotModule::setSubtitle(const char *subtitle)
 }
 
 
+void
+AbstractPlotModule::setSubtitle(const std::string &subtitle)
+{
+    impl_->subtitle_ = subtitle;
+}
+
+
 void
 AbstractPlotModule::setXLabel(const char *label)
 {
index 3bd46b2a51c2ec8f36907ce53c7f741bea79cc08..78c48d5eecfe2957b58de860703487e8b7ff9e4f 100644 (file)
@@ -177,10 +177,14 @@ class AbstractPlotModule : public AnalysisDataModuleSerial
          * Set plot title.
          */
         void setTitle(const char *title);
+        //! \copydoc setTitle(const char *)
+        void setTitle(const std::string &title);
         /*! \brief
          * Set plot subtitle.
          */
         void setSubtitle(const char *subtitle);
+        //! \copydoc setSubtitle(const char *)
+        void setSubtitle(const std::string &subtitle);
         /*! \brief
          * Set X axis label.
          */
index acac0e9bebee62e18c135cd139f45f68d0237ad0..46760d8f19dfa1ae3a86a6cb540328cc9a50a40e 100644 (file)
@@ -189,9 +189,6 @@ gmx_principal(int argc, char *argv[]);
 int
 gmx_rama(int argc, char *argv[]);
 
-int
-gmx_rdf(int argc, char *argv[]);
-
 int
 gmx_rotmat(int argc, char *argv[]);
 
diff --git a/src/gromacs/gmxana/gmx_rdf.c b/src/gromacs/gmxana/gmx_rdf.c
deleted file mode 100644 (file)
index 067483b..0000000
+++ /dev/null
@@ -1,946 +0,0 @@
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,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.
- */
-
-#include "gmxpre.h"
-
-#include <math.h>
-
-#include "gromacs/commandline/pargs.h"
-#include "gromacs/correlationfunctions/integrate.h"
-#include "gromacs/fileio/tpxio.h"
-#include "gromacs/fileio/trxio.h"
-#include "gromacs/fileio/xvgr.h"
-#include "gromacs/gmxana/gmx_ana.h"
-#include "gromacs/gmxana/gstat.h"
-#include "gromacs/legacyheaders/macros.h"
-#include "gromacs/legacyheaders/names.h"
-#include "gromacs/legacyheaders/nrnb.h"
-#include "gromacs/legacyheaders/typedefs.h"
-#include "gromacs/legacyheaders/viewit.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/pbcutil/pbc.h"
-#include "gromacs/pbcutil/rmpbc.h"
-#include "gromacs/topology/index.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/futil.h"
-#include "gromacs/utility/smalloc.h"
-
-static void check_box_c(matrix box)
-{
-    if (fabs(box[ZZ][XX]) > GMX_REAL_EPS*box[ZZ][ZZ] ||
-        fabs(box[ZZ][YY]) > GMX_REAL_EPS*box[ZZ][ZZ])
-    {
-        gmx_fatal(FARGS,
-                  "The last box vector is not parallel to the z-axis: %f %f %f",
-                  box[ZZ][XX], box[ZZ][YY], box[ZZ][ZZ]);
-    }
-}
-
-static void calc_comg(int is, int *coi, int *index, gmx_bool bMass, t_atom *atom,
-                      rvec *x, rvec *x_comg)
-{
-    int  c, i, d;
-    rvec xc;
-    real mtot, m;
-
-    if (bMass && atom == NULL)
-    {
-        gmx_fatal(FARGS, "No masses available while mass weighting was requested");
-    }
-
-    for (c = 0; c < is; c++)
-    {
-        clear_rvec(xc);
-        mtot = 0;
-        for (i = coi[c]; i < coi[c+1]; i++)
-        {
-            if (bMass)
-            {
-                m = atom[index[i]].m;
-                for (d = 0; d < DIM; d++)
-                {
-                    xc[d] += m*x[index[i]][d];
-                }
-                mtot += m;
-            }
-            else
-            {
-                rvec_inc(xc, x[index[i]]);
-                mtot += 1.0;
-            }
-        }
-        svmul(1/mtot, xc, x_comg[c]);
-    }
-}
-
-static void split_group(int isize, int *index, char *grpname,
-                        t_topology *top, char type,
-                        int *is_out, int **coi_out)
-{
-    t_block *mols = NULL;
-    t_atom  *atom = NULL;
-    int      is, *coi;
-    int      cur, mol, res, i, a, i1;
-
-    /* Split up the group in molecules or residues */
-    switch (type)
-    {
-        case 'm':
-            mols = &top->mols;
-            break;
-        case 'r':
-            atom = top->atoms.atom;
-            break;
-        default:
-            gmx_fatal(FARGS, "Unknown rdf option '%s'", type);
-    }
-    snew(coi, isize+1);
-    is  = 0;
-    cur = -1;
-    mol = 0;
-    for (i = 0; i < isize; i++)
-    {
-        a = index[i];
-        if (type == 'm')
-        {
-            /* Check if the molecule number has changed */
-            i1 = mols->index[mol+1];
-            while (a >= i1)
-            {
-                mol++;
-                i1 = mols->index[mol+1];
-            }
-            if (mol != cur)
-            {
-                coi[is++] = i;
-                cur       = mol;
-            }
-        }
-        else if (type == 'r')
-        {
-            /* Check if the residue index has changed */
-            res = atom[a].resind;
-            if (res != cur)
-            {
-                coi[is++] = i;
-                cur       = res;
-            }
-        }
-    }
-    coi[is] = i;
-    srenew(coi, is+1);
-    printf("Group '%s' of %d atoms consists of %d %s\n",
-           grpname, isize, is,
-           (type == 'm' ? "molecules" : "residues"));
-
-    *is_out  = is;
-    *coi_out = coi;
-}
-
-static void do_rdf(const char *fnNDX, const char *fnTPS, const char *fnTRX,
-                   const char *fnRDF, const char *fnCNRDF, const char *fnHQ,
-                   gmx_bool bCM, const char *close,
-                   const char **rdft, gmx_bool bXY, gmx_bool bPBC, gmx_bool bNormalize,
-                   real cutoff, real binwidth, real fade, int ng,
-                   const output_env_t oenv)
-{
-    FILE          *fp;
-    t_trxstatus   *status;
-    char           outf1[STRLEN], outf2[STRLEN];
-    char           title[STRLEN], gtitle[STRLEN], refgt[30];
-    int            g, natoms, i, ii, j, k, nbin, j0, j1, n, nframes;
-    int          **count;
-    char         **grpname;
-    int           *isize, isize_cm = 0, nrdf = 0, max_i, isize0, isize_g;
-    atom_id      **index, *index_cm = NULL;
-    gmx_int64_t   *sum;
-    real           t, rmax2, cut2, r, r2, r2ii, invhbinw, normfac;
-    real           segvol, spherevol, prev_spherevol, **rdf;
-    rvec          *x, dx, *x0 = NULL, *x_i1, xi;
-    real          *inv_segvol, invvol, invvol_sum, rho;
-    gmx_bool       bClose, *bExcl, bTop, bNonSelfExcl;
-    matrix         box, box_pbc;
-    int          **npairs;
-    atom_id        ix, jx, ***pairs;
-    t_topology    *top  = NULL;
-    int            ePBC = -1, ePBCrdf = -1;
-    t_block       *mols = NULL;
-    t_blocka      *excl;
-    t_atom        *atom = NULL;
-    t_pbc          pbc;
-    gmx_rmpbc_t    gpbc = NULL;
-    int           *is   = NULL, **coi = NULL, cur, mol, i1, res, a;
-
-    excl = NULL;
-
-    bClose = (close[0] != 'n');
-
-    if (fnTPS)
-    {
-        snew(top, 1);
-        bTop = read_tps_conf(fnTPS, title, top, &ePBC, &x, NULL, box, TRUE);
-        if (bTop && !(bCM || bClose))
-        {
-            /* get exclusions from topology */
-            excl = &(top->excls);
-        }
-    }
-    snew(grpname, ng+1);
-    snew(isize, ng+1);
-    snew(index, ng+1);
-    fprintf(stderr, "\nSelect a reference group and %d group%s\n",
-            ng, ng == 1 ? "" : "s");
-    if (fnTPS)
-    {
-        get_index(&(top->atoms), fnNDX, ng+1, isize, index, grpname);
-        atom = top->atoms.atom;
-    }
-    else
-    {
-        rd_index(fnNDX, ng+1, isize, index, grpname);
-    }
-
-    if (bCM || bClose)
-    {
-        snew(is, 1);
-        snew(coi, 1);
-        if (bClose)
-        {
-            split_group(isize[0], index[0], grpname[0], top, close[0], &is[0], &coi[0]);
-        }
-    }
-    if (rdft[0][0] != 'a')
-    {
-        /* Split up all the groups in molecules or residues */
-        srenew(is, ng+1);
-        srenew(coi, ng+1);
-        for (g = ((bCM || bClose) ? 1 : 0); g < ng+1; g++)
-        {
-            split_group(isize[g], index[g], grpname[g], top, rdft[0][0], &is[g], &coi[g]);
-        }
-    }
-
-    if (bCM)
-    {
-        is[0] = 1;
-        snew(coi[0], is[0]+1);
-        coi[0][0] = 0;
-        coi[0][1] = isize[0];
-        isize0    = is[0];
-        snew(x0, isize0);
-    }
-    else if (bClose || rdft[0][0] != 'a')
-    {
-        isize0 = is[0];
-        snew(x0, isize0);
-    }
-    else
-    {
-        isize0 = isize[0];
-    }
-
-    natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
-    if (!natoms)
-    {
-        gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
-    }
-    if (fnTPS)
-    {
-        /* check with topology */
-        if (natoms > top->atoms.nr)
-        {
-            gmx_fatal(FARGS, "Trajectory (%d atoms) does not match topology (%d atoms)",
-                      natoms, top->atoms.nr);
-        }
-    }
-    /* check with index groups */
-    for (i = 0; i < ng+1; i++)
-    {
-        for (j = 0; j < isize[i]; j++)
-        {
-            if (index[i][j] >= natoms)
-            {
-                gmx_fatal(FARGS, "Atom index (%d) in index group %s (%d atoms) larger "
-                          "than number of atoms in trajectory (%d atoms)",
-                          index[i][j], grpname[i], isize[i], natoms);
-            }
-        }
-    }
-
-    /* initialize some handy things */
-    if (ePBC == -1)
-    {
-        ePBC = guess_ePBC(box);
-    }
-    copy_mat(box, box_pbc);
-    if (bXY)
-    {
-        check_box_c(box);
-        switch (ePBC)
-        {
-            case epbcXYZ:
-            case epbcXY:   ePBCrdf = epbcXY;   break;
-            case epbcNONE: ePBCrdf = epbcNONE; break;
-            default:
-                gmx_fatal(FARGS, "xy-rdf's are not supported for pbc type'%s'",
-                          EPBC(ePBC));
-                break;
-        }
-        /* Make sure the z-height does not influence the cut-off */
-        box_pbc[ZZ][ZZ] = 2*max(box[XX][XX], box[YY][YY]);
-    }
-    else
-    {
-        ePBCrdf = ePBC;
-    }
-    if (bPBC)
-    {
-        rmax2   = 0.99*0.99*max_cutoff2(bXY ? epbcXY : epbcXYZ, box_pbc);
-    }
-    else
-    {
-        rmax2   = sqr(3*max(box[XX][XX], max(box[YY][YY], box[ZZ][ZZ])));
-    }
-    if (debug)
-    {
-        fprintf(debug, "rmax2 = %g\n", rmax2);
-    }
-
-    /* We use the double amount of bins, so we can correctly
-     * write the rdf and rdf_cn output at i*binwidth values.
-     */
-    nbin     = (int)(sqrt(rmax2) * 2 / binwidth);
-    invhbinw = 2.0 / binwidth;
-    cut2     = sqr(cutoff);
-
-    snew(count, ng);
-    snew(pairs, ng);
-    snew(npairs, ng);
-
-    snew(bExcl, natoms);
-    max_i = 0;
-    for (g = 0; g < ng; g++)
-    {
-        if (isize[g+1] > max_i)
-        {
-            max_i = isize[g+1];
-        }
-
-        /* this is THE array */
-        snew(count[g], nbin+1);
-
-        /* make pairlist array for groups and exclusions */
-        snew(pairs[g], isize[0]);
-        snew(npairs[g], isize[0]);
-        for (i = 0; i < isize[0]; i++)
-        {
-            /* We can only have exclusions with atomic rdfs */
-            if (!(bCM || bClose || rdft[0][0] != 'a'))
-            {
-                ix = index[0][i];
-                for (j = 0; j < natoms; j++)
-                {
-                    bExcl[j] = FALSE;
-                }
-                /* exclusions? */
-                if (excl)
-                {
-                    for (j = excl->index[ix]; j < excl->index[ix+1]; j++)
-                    {
-                        bExcl[excl->a[j]] = TRUE;
-                    }
-                }
-                k = 0;
-                snew(pairs[g][i], isize[g+1]);
-                bNonSelfExcl = FALSE;
-                for (j = 0; j < isize[g+1]; j++)
-                {
-                    jx = index[g+1][j];
-                    if (!bExcl[jx])
-                    {
-                        pairs[g][i][k++] = jx;
-                    }
-                    else if (ix != jx)
-                    {
-                        /* Check if we have exclusions other than self exclusions */
-                        bNonSelfExcl = TRUE;
-                    }
-                }
-                if (bNonSelfExcl)
-                {
-                    npairs[g][i] = k;
-                    srenew(pairs[g][i], npairs[g][i]);
-                }
-                else
-                {
-                    /* Save a LOT of memory and some cpu cycles */
-                    npairs[g][i] = -1;
-                    sfree(pairs[g][i]);
-                }
-            }
-            else
-            {
-                npairs[g][i] = -1;
-            }
-        }
-    }
-    sfree(bExcl);
-
-    snew(x_i1, max_i);
-    nframes    = 0;
-    invvol_sum = 0;
-    if (bPBC && (NULL != top))
-    {
-        gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
-    }
-
-    do
-    {
-        /* Must init pbc every step because of pressure coupling */
-        copy_mat(box, box_pbc);
-        if (bPBC)
-        {
-            if (top != NULL)
-            {
-                gmx_rmpbc(gpbc, natoms, box, x);
-            }
-            if (bXY)
-            {
-                check_box_c(box);
-                clear_rvec(box_pbc[ZZ]);
-            }
-            set_pbc(&pbc, ePBCrdf, box_pbc);
-
-            if (bXY)
-            {
-                /* Set z-size to 1 so we get the surface iso the volume */
-                box_pbc[ZZ][ZZ] = 1;
-            }
-        }
-        invvol      = 1/det(box_pbc);
-        invvol_sum += invvol;
-
-        if (bCM)
-        {
-            /* Calculate center of mass of the whole group */
-            calc_comg(is[0], coi[0], index[0], TRUE, atom, x, x0);
-        }
-        else if (!bClose && rdft[0][0] != 'a')
-        {
-            calc_comg(is[0], coi[0], index[0], rdft[0][6] == 'm', atom, x, x0);
-        }
-
-        for (g = 0; g < ng; g++)
-        {
-            if (rdft[0][0] == 'a')
-            {
-                /* Copy the indexed coordinates to a continuous array */
-                for (i = 0; i < isize[g+1]; i++)
-                {
-                    copy_rvec(x[index[g+1][i]], x_i1[i]);
-                }
-            }
-            else
-            {
-                /* Calculate the COMs/COGs and store in x_i1 */
-                calc_comg(is[g+1], coi[g+1], index[g+1], rdft[0][6] == 'm', atom, x, x_i1);
-            }
-
-            for (i = 0; i < isize0; i++)
-            {
-                if (bClose)
-                {
-                    /* Special loop, since we need to determine the minimum distance
-                     * over all selected atoms in the reference molecule/residue.
-                     */
-                    if (rdft[0][0] == 'a')
-                    {
-                        isize_g = isize[g+1];
-                    }
-                    else
-                    {
-                        isize_g = is[g+1];
-                    }
-                    for (j = 0; j < isize_g; j++)
-                    {
-                        r2 = 1e30;
-                        /* Loop over the selected atoms in the reference molecule */
-                        for (ii = coi[0][i]; ii < coi[0][i+1]; ii++)
-                        {
-                            if (bPBC)
-                            {
-                                pbc_dx(&pbc, x[index[0][ii]], x_i1[j], dx);
-                            }
-                            else
-                            {
-                                rvec_sub(x[index[0][ii]], x_i1[j], dx);
-                            }
-                            if (bXY)
-                            {
-                                r2ii = dx[XX]*dx[XX] + dx[YY]*dx[YY];
-                            }
-                            else
-                            {
-                                r2ii = iprod(dx, dx);
-                            }
-                            if (r2ii < r2)
-                            {
-                                r2 = r2ii;
-                            }
-                        }
-                        if (r2 > cut2 && r2 <= rmax2)
-                        {
-                            count[g][(int)(sqrt(r2)*invhbinw)]++;
-                        }
-                    }
-                }
-                else
-                {
-                    /* Real rdf between points in space */
-                    if (bCM || rdft[0][0] != 'a')
-                    {
-                        copy_rvec(x0[i], xi);
-                    }
-                    else
-                    {
-                        copy_rvec(x[index[0][i]], xi);
-                    }
-                    if (rdft[0][0] == 'a' && npairs[g][i] >= 0)
-                    {
-                        /* Expensive loop, because of indexing */
-                        for (j = 0; j < npairs[g][i]; j++)
-                        {
-                            jx = pairs[g][i][j];
-                            if (bPBC)
-                            {
-                                pbc_dx(&pbc, xi, x[jx], dx);
-                            }
-                            else
-                            {
-                                rvec_sub(xi, x[jx], dx);
-                            }
-
-                            if (bXY)
-                            {
-                                r2 = dx[XX]*dx[XX] + dx[YY]*dx[YY];
-                            }
-                            else
-                            {
-                                r2 = iprod(dx, dx);
-                            }
-                            if (r2 > cut2 && r2 <= rmax2)
-                            {
-                                count[g][(int)(sqrt(r2)*invhbinw)]++;
-                            }
-                        }
-                    }
-                    else
-                    {
-                        /* Cheaper loop, no exclusions */
-                        if (rdft[0][0] == 'a')
-                        {
-                            isize_g = isize[g+1];
-                        }
-                        else
-                        {
-                            isize_g = is[g+1];
-                        }
-                        for (j = 0; j < isize_g; j++)
-                        {
-                            if (bPBC)
-                            {
-                                pbc_dx(&pbc, xi, x_i1[j], dx);
-                            }
-                            else
-                            {
-                                rvec_sub(xi, x_i1[j], dx);
-                            }
-                            if (bXY)
-                            {
-                                r2 = dx[XX]*dx[XX] + dx[YY]*dx[YY];
-                            }
-                            else
-                            {
-                                r2 = iprod(dx, dx);
-                            }
-                            if (r2 > cut2 && r2 <= rmax2)
-                            {
-                                count[g][(int)(sqrt(r2)*invhbinw)]++;
-                            }
-                        }
-                    }
-                }
-            }
-        }
-        nframes++;
-    }
-    while (read_next_x(oenv, status, &t, x, box));
-    fprintf(stderr, "\n");
-
-    if (bPBC && (NULL != top))
-    {
-        gmx_rmpbc_done(gpbc);
-    }
-
-    close_trj(status);
-
-    sfree(x);
-
-    /* Average volume */
-    invvol = invvol_sum/nframes;
-
-    /* Calculate volume of sphere segments or length of circle segments */
-    snew(inv_segvol, (nbin+1)/2);
-    prev_spherevol = 0;
-    for (i = 0; (i < (nbin+1)/2); i++)
-    {
-        r = (i + 0.5)*binwidth;
-        if (bXY)
-        {
-            spherevol = M_PI*r*r;
-        }
-        else
-        {
-            spherevol = (4.0/3.0)*M_PI*r*r*r;
-        }
-        segvol         = spherevol-prev_spherevol;
-        inv_segvol[i]  = 1.0/segvol;
-        prev_spherevol = spherevol;
-    }
-
-    snew(rdf, ng);
-    for (g = 0; g < ng; g++)
-    {
-        /* We have to normalize by dividing by the number of frames */
-        if (rdft[0][0] == 'a')
-        {
-            normfac = 1.0/(nframes*invvol*isize0*isize[g+1]);
-        }
-        else
-        {
-            normfac = 1.0/(nframes*invvol*isize0*is[g+1]);
-        }
-
-        /* Do the normalization */
-        nrdf = max((nbin+1)/2, 1+2*fade/binwidth);
-        snew(rdf[g], nrdf);
-        for (i = 0; i < (nbin+1)/2; i++)
-        {
-            r = i*binwidth;
-            if (i == 0)
-            {
-                j = count[g][0];
-            }
-            else
-            {
-                j = count[g][i*2-1] + count[g][i*2];
-            }
-            if ((fade > 0) && (r >= fade))
-            {
-                rdf[g][i] = 1 + (j*inv_segvol[i]*normfac-1)*exp(-16*sqr(r/fade-1));
-            }
-            else
-            {
-                if (bNormalize)
-                {
-                    rdf[g][i] = j*inv_segvol[i]*normfac;
-                }
-                else
-                {
-                    rdf[g][i] = j/(binwidth*isize0*nframes);
-                }
-            }
-        }
-        for (; (i < nrdf); i++)
-        {
-            rdf[g][i] = 1.0;
-        }
-    }
-
-    if (rdft[0][0] == 'a')
-    {
-        sprintf(gtitle, "Radial distribution");
-    }
-    else
-    {
-        sprintf(gtitle, "Radial distribution of %s %s",
-                rdft[0][0] == 'm' ? "molecule" : "residue",
-                rdft[0][6] == 'm' ? "COM" : "COG");
-    }
-    fp = xvgropen(fnRDF, gtitle, "r", "", oenv);
-    if (bCM)
-    {
-        sprintf(refgt, " %s", "COM");
-    }
-    else if (bClose)
-    {
-        sprintf(refgt, " closest atom in %s.", close);
-    }
-    else
-    {
-        sprintf(refgt, "%s", "");
-    }
-    if (ng == 1)
-    {
-        if (output_env_get_print_xvgr_codes(oenv))
-        {
-            fprintf(fp, "@ subtitle \"%s%s - %s\"\n", grpname[0], refgt, grpname[1]);
-        }
-    }
-    else
-    {
-        if (output_env_get_print_xvgr_codes(oenv))
-        {
-            fprintf(fp, "@ subtitle \"reference %s%s\"\n", grpname[0], refgt);
-        }
-        xvgr_legend(fp, ng, (const char**)(grpname+1), oenv);
-    }
-    for (i = 0; (i < nrdf); i++)
-    {
-        fprintf(fp, "%10g", i*binwidth);
-        for (g = 0; g < ng; g++)
-        {
-            fprintf(fp, " %10g", rdf[g][i]);
-        }
-        fprintf(fp, "\n");
-    }
-    gmx_ffclose(fp);
-
-    do_view(oenv, fnRDF, NULL);
-
-    /* h(Q) function: fourier transform of rdf */
-    if (fnHQ)
-    {
-        int   nhq = 401;
-        real *hq, *integrand, Q;
-
-        /* Get a better number density later! */
-        rho = isize[1]*invvol;
-        snew(hq, nhq);
-        snew(integrand, nrdf);
-        for (i = 0; (i < nhq); i++)
-        {
-            Q            = i*0.5;
-            integrand[0] = 0;
-            for (j = 1; (j < nrdf); j++)
-            {
-                r             = j*binwidth;
-                integrand[j]  = (Q == 0) ? 1.0 : sin(Q*r)/(Q*r);
-                integrand[j] *= 4.0*M_PI*rho*r*r*(rdf[0][j]-1.0);
-            }
-            hq[i] = print_and_integrate(debug, nrdf, binwidth, integrand, NULL, 0);
-        }
-        fp = xvgropen(fnHQ, "h(Q)", "Q(/nm)", "h(Q)", oenv);
-        for (i = 0; (i < nhq); i++)
-        {
-            fprintf(fp, "%10g %10g\n", i*0.5, hq[i]);
-        }
-        gmx_ffclose(fp);
-        do_view(oenv, fnHQ, NULL);
-        sfree(hq);
-        sfree(integrand);
-    }
-
-    if (fnCNRDF)
-    {
-        normfac = 1.0/(isize0*nframes);
-        fp      = xvgropen(fnCNRDF, "Cumulative Number RDF", "r", "number", oenv);
-        if (ng == 1)
-        {
-            if (output_env_get_print_xvgr_codes(oenv))
-            {
-                fprintf(fp, "@ subtitle \"%s-%s\"\n", grpname[0], grpname[1]);
-            }
-        }
-        else
-        {
-            if (output_env_get_print_xvgr_codes(oenv))
-            {
-                fprintf(fp, "@ subtitle \"reference %s\"\n", grpname[0]);
-            }
-            xvgr_legend(fp, ng, (const char**)(grpname+1), oenv);
-        }
-        snew(sum, ng);
-        for (i = 0; (i <= nbin/2); i++)
-        {
-            fprintf(fp, "%10g", i*binwidth);
-            for (g = 0; g < ng; g++)
-            {
-                fprintf(fp, " %10g", (real)((double)sum[g]*normfac));
-                if (i*2+1 < nbin)
-                {
-                    sum[g] += count[g][i*2] + count[g][i*2+1];
-                }
-            }
-            fprintf(fp, "\n");
-        }
-        gmx_ffclose(fp);
-        sfree(sum);
-
-        do_view(oenv, fnCNRDF, NULL);
-    }
-
-    for (g = 0; g < ng; g++)
-    {
-        sfree(rdf[g]);
-    }
-    sfree(rdf);
-}
-
-
-int gmx_rdf(int argc, char *argv[])
-{
-    const char        *desc[] = {
-        "The structure of liquids can be studied by either neutron or X-ray",
-        "scattering. The most common way to describe liquid structure is by a",
-        "radial distribution function. However, this is not easy to obtain from",
-        "a scattering experiment.[PAR]",
-        "[THISMODULE] calculates radial distribution functions in different ways.",
-        "The normal method is around a (set of) particle(s), the other methods",
-        "are around the center of mass of a set of particles ([TT]-com[tt])",
-        "or to the closest particle in a set ([TT]-surf[tt]).",
-        "With all methods, the RDF can also be calculated around axes parallel",
-        "to the [IT]z[it]-axis with option [TT]-xy[tt].",
-        "With option [TT]-surf[tt] normalization can not be used.[PAR]",
-        "The option [TT]-rdf[tt] sets the type of RDF to be computed.",
-        "Default is for atoms or particles, but one can also select center",
-        "of mass or geometry of molecules or residues. In all cases, only",
-        "the atoms in the index groups are taken into account.",
-        "For molecules and/or the center of mass option, a run input file",
-        "is required.",
-        "Weighting other than COM or COG can currently only be achieved",
-        "by providing a run input file with different masses.",
-        "Options [TT]-com[tt] and [TT]-surf[tt] also work in conjunction",
-        "with [TT]-rdf[tt].[PAR]",
-        "If a run input file is supplied ([TT]-s[tt]) and [TT]-rdf[tt] is set",
-        "to [TT]atom[tt], exclusions defined",
-        "in that file are taken into account when calculating the RDF.",
-        "The option [TT]-cut[tt] is meant as an alternative way to avoid",
-        "intramolecular peaks in the RDF plot.",
-        "It is however better to supply a run input file with a higher number of",
-        "exclusions. For e.g. benzene a topology, setting nrexcl to 5",
-        "would eliminate all intramolecular contributions to the RDF.",
-        "Note that all atoms in the selected groups are used, also the ones",
-        "that don't have Lennard-Jones interactions.[PAR]",
-        "Option [TT]-cn[tt] produces the cumulative number RDF,",
-        "i.e. the average number of particles within a distance r.[PAR]"
-    };
-    static gmx_bool    bCM     = FALSE, bXY = FALSE, bPBC = TRUE, bNormalize = TRUE;
-    static real        cutoff  = 0, binwidth = 0.002, fade = 0.0;
-    static int         ngroups = 1;
-
-    static const char *closet[] = { NULL, "no", "mol", "res", NULL };
-    static const char *rdft[]   = { NULL, "atom", "mol_com", "mol_cog", "res_com", "res_cog", NULL };
-
-    t_pargs            pa[] = {
-        { "-bin",      FALSE, etREAL, {&binwidth},
-          "Binwidth (nm)" },
-        { "-com",      FALSE, etBOOL, {&bCM},
-          "RDF with respect to the center of mass of first group" },
-        { "-surf",     FALSE, etENUM, {closet},
-          "RDF with respect to the surface of the first group" },
-        { "-rdf",   FALSE, etENUM, {rdft},
-          "RDF type" },
-        { "-pbc",      FALSE, etBOOL, {&bPBC},
-          "Use periodic boundary conditions for computing distances. Without PBC the maximum range will be three times the largest box edge." },
-        { "-norm",     FALSE, etBOOL, {&bNormalize},
-          "Normalize for volume and density" },
-        { "-xy",       FALSE, etBOOL, {&bXY},
-          "Use only the x and y components of the distance" },
-        { "-cut",      FALSE, etREAL, {&cutoff},
-          "Shortest distance (nm) to be considered"},
-        { "-ng",       FALSE, etINT, {&ngroups},
-          "Number of secondary groups to compute RDFs around a central group" },
-        { "-fade",     FALSE, etREAL, {&fade},
-          "From this distance onwards the RDF is tranformed by g'(r) = 1 + [g(r)-1] exp(-(r/fade-1)^2 to make it go to 1 smoothly. If fade is 0.0 nothing is done." }
-    };
-#define NPA asize(pa)
-    const char        *fnTPS, *fnNDX;
-    output_env_t       oenv;
-
-    t_filenm           fnm[] = {
-        { efTRX, "-f",  NULL,     ffREAD },
-        { efTPS, NULL,  NULL,     ffOPTRD },
-        { efNDX, NULL,  NULL,     ffOPTRD },
-        { efXVG, "-o",  "rdf",    ffWRITE },
-        { efXVG, "-cn", "rdf_cn", ffOPTWR },
-        { efXVG, "-hq", "hq",     ffOPTWR },
-    };
-#define NFILE asize(fnm)
-    if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
-                           NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
-    {
-        return 0;
-    }
-
-    if (bCM || closet[0][0] != 'n' || rdft[0][0] == 'm' || rdft[0][6] == 'm')
-    {
-        fnTPS = ftp2fn(efTPS, NFILE, fnm);
-    }
-    else
-    {
-        fnTPS = ftp2fn_null(efTPS, NFILE, fnm);
-    }
-    fnNDX = ftp2fn_null(efNDX, NFILE, fnm);
-
-    if (!fnTPS && !fnNDX)
-    {
-        gmx_fatal(FARGS, "Neither index file nor topology file specified\n"
-                  "Nothing to do!");
-    }
-
-    if (closet[0][0] != 'n')
-    {
-        if (bCM)
-        {
-            gmx_fatal(FARGS, "Can not have both -com and -surf");
-        }
-        if (bNormalize)
-        {
-            fprintf(stderr, "Turning of normalization because of option -surf\n");
-            bNormalize = FALSE;
-        }
-    }
-
-    do_rdf(fnNDX, fnTPS, ftp2fn(efTRX, NFILE, fnm),
-           opt2fn("-o", NFILE, fnm), opt2fn_null("-cn", NFILE, fnm),
-           opt2fn_null("-hq", NFILE, fnm),
-           bCM, closet[0], rdft, bXY, bPBC, bNormalize, cutoff, binwidth, fade, ngroups,
-           oenv);
-
-    return 0;
-}
index df6c99f23bc7345e9dcef00ec07d54c5828926cc..7fc2601f6df577ae02d2d9245c960d2316655d96 100644 (file)
@@ -50,6 +50,7 @@
 #include "modules/distance.h"
 #include "modules/freevolume.h"
 #include "modules/pairdist.h"
+#include "modules/rdf.h"
 #include "modules/sasa.h"
 #include "modules/select.h"
 
@@ -92,6 +93,7 @@ void registerTrajectoryAnalysisModules(CommandLineModuleManager *manager)
     registerModule<DistanceInfo>(manager, group);
     registerModule<FreeVolumeInfo>(manager, group);
     registerModule<PairDistanceInfo>(manager, group);
+    registerModule<RdfInfo>(manager, group);
     registerModule<SasaInfo>(manager, group);
     registerModule<SelectInfo>(manager, group);
 }
diff --git a/src/gromacs/trajectoryanalysis/modules/rdf.cpp b/src/gromacs/trajectoryanalysis/modules/rdf.cpp
new file mode 100644 (file)
index 0000000..468b729
--- /dev/null
@@ -0,0 +1,663 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,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
+ * Implements gmx::analysismodules::Rdf.
+ *
+ * \author Teemu Murtola <teemu.murtola@gmail.com> (C++ conversion)
+ * \ingroup module_trajectoryanalysis
+ */
+#include "gmxpre.h"
+
+#include "rdf.h"
+
+#include <cmath>
+
+#include <algorithm>
+#include <limits>
+#include <string>
+#include <vector>
+
+#include "gromacs/analysisdata/analysisdata.h"
+#include "gromacs/analysisdata/modules/average.h"
+#include "gromacs/analysisdata/modules/histogram.h"
+#include "gromacs/analysisdata/modules/plot.h"
+#include "gromacs/fileio/trx.h"
+#include "gromacs/math/utilities.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/options/basicoptions.h"
+#include "gromacs/options/filenameoption.h"
+#include "gromacs/options/options.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/selection/nbsearch.h"
+#include "gromacs/selection/selection.h"
+#include "gromacs/selection/selectionoption.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/trajectoryanalysis/analysismodule.h"
+#include "gromacs/trajectoryanalysis/analysissettings.h"
+#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/stringutil.h"
+
+namespace gmx
+{
+
+namespace analysismodules
+{
+
+namespace
+{
+
+//! \addtogroup module_trajectoryanalysis
+//! \{
+
+/********************************************************************
+ * Actual analysis module
+ */
+
+/*! \brief
+ * Implements `gmx rdf` trajectory analysis module.
+ */
+class Rdf : public TrajectoryAnalysisModule
+{
+    public:
+        Rdf();
+
+        virtual void initOptions(Options                    *options,
+                                 TrajectoryAnalysisSettings *settings);
+        virtual void optionsFinished(Options                    *options,
+                                     TrajectoryAnalysisSettings *settings);
+        virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
+                                  const TopologyInformation        &top);
+        virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
+                                         const t_trxframe                 &fr);
+
+        virtual TrajectoryAnalysisModuleDataPointer startFrames(
+            const AnalysisDataParallelOptions &opt,
+            const SelectionCollection         &selections);
+        virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
+                                  TrajectoryAnalysisModuleData *pdata);
+
+        virtual void finishAnalysis(int nframes);
+        virtual void writeOutput();
+
+    private:
+        std::string                               fnRdf_;
+        std::string                               fnCumulative_;
+        std::string                               surface_;
+        AnalysisDataPlotSettings                  plotSettings_;
+
+        /*! \brief
+         * Reference selection to compute RDFs around.
+         *
+         * With -surf, Selection::originalIds() and Selection::mappedIds()
+         * store the index of the surface group to which that position belongs.
+         * The RDF is computed by finding the nearest position from each
+         * surface group for each position, and then binning those distances.
+         */
+        Selection                                 refSel_;
+        /*! \brief
+         * Selections to compute RDFs for.
+         */
+        SelectionList                             sel_;
+
+        /*! \brief
+         * Raw pairwise distance data from which the RDF is computed.
+         *
+         * There is a data set for each selection in `sel_`, with a single
+         * column.  Each point set will contain a single pairwise distance
+         * that contributes to the RDF.
+         */
+        AnalysisData                              pairDist_;
+        /*! \brief
+         * Normalization factors for each frame.
+         *
+         * The first column contains the number of positions in `refSel_` for
+         * that frame (with surface RDF, the number of groups).  There are
+         * `sel_.size()` more columns, each containing the number density of
+         * positions for one selection.
+         */
+        AnalysisData                              normFactors_;
+        /*! \brief
+         * Histogram module that computes the actual RDF from `pairDist_`.
+         *
+         * The per-frame histograms are raw pair counts in each bin;
+         * the averager is normalized by the average number of reference
+         * positions (average of the first column of `normFactors_`).
+         */
+        AnalysisDataSimpleHistogramModulePointer  pairCounts_;
+        /*! \brief
+         * Average normalization factors.
+         */
+        AnalysisDataAverageModulePointer          normAve_;
+        //! Neighborhood search with `refSel_` as the reference positions.
+        AnalysisNeighborhood                      nb_;
+
+        // User input options.
+        double                                    binwidth_;
+        double                                    cutoff_;
+        double                                    rmax_;
+        bool                                      bNormalize_;
+        bool                                      bXY_;
+        bool                                      bExclusions_;
+
+        // Pre-computed values for faster access during analysis.
+        real                                      cut2_;
+        real                                      rmax2_;
+        int                                       surfaceGroupCount_;
+
+        // Copy and assign disallowed by base.
+};
+
+Rdf::Rdf()
+    : TrajectoryAnalysisModule(RdfInfo::name, RdfInfo::shortDescription),
+      pairCounts_(new AnalysisDataSimpleHistogramModule()),
+      normAve_(new AnalysisDataAverageModule()),
+      binwidth_(0.002), cutoff_(0.0), rmax_(0.0),
+      bNormalize_(true), bXY_(false), bExclusions_(false),
+      cut2_(0.0), rmax2_(0.0), surfaceGroupCount_(0)
+{
+    pairDist_.setMultipoint(true);
+    pairDist_.addModule(pairCounts_);
+    registerAnalysisDataset(&pairDist_, "pairdist");
+    registerBasicDataset(pairCounts_.get(), "paircount");
+
+    normFactors_.addModule(normAve_);
+    registerAnalysisDataset(&normFactors_, "norm");
+}
+
+void
+Rdf::initOptions(Options *options, TrajectoryAnalysisSettings * /*settings*/)
+{
+    static const char *const desc[] = {
+        "[THISMODULE] calculates radial distribution functions from one",
+        "refernce set of position (set with [TT]-ref[tt]) to one or more",
+        "sets of positions (set with [TT]-sel[tt]).  To compute the RDF with",
+        "respect to the closest position in a set in [TT]-ref[tt] instead, use",
+        "[TT]-surf[tt]: if set, then [TT]-ref[tt] is partitioned into sets",
+        "based on the value of [TT]-surf[tt], and the closest position in each",
+        "set is used. To compute the RDF around axes parallel to the",
+        "[IT]z[it]-axis, i.e., only in the [IT]x[it]-[IT]y[it] plane, use",
+        "[TT]-xy[tt].[PAR]",
+        "To set the bin width and maximum distance to use in the RDF, use",
+        "[TT]-bin[tt] and [TT]-rmax[tt], respectively. The latter can be",
+        "used to limit the computational cost if the RDF is not of interest",
+        "up to the default (half of the box size with PBC, three times the",
+        "box size without PBC).[PAR]",
+        "To use exclusions from the topology ([TT]-s[tt]), set [TT]-excl[tt]",
+        "and ensure that both [TT]-ref[tt] and [TT]-sel[tt] only select atoms.",
+        "A rougher alternative to exclude intra-molecular peaks is to set",
+        "[TT]-cut[tt] to a non-zero value to clear the RDF at small",
+        "distances.[PAR]",
+        "The RDFs are normalized by 1) average number of positions in",
+        "[TT]-ref[tt] (the number of groups with [TT]-surf[tt]), 2) volume",
+        "of the bin, and 3) average particle density of [TT]-sel[tt] positions",
+        "for that selection. To only use the first factor for normalization,",
+        "set [TT]-nonorm[tt]. In this case, the RDF is only scaled with the",
+        "bin width to make the integral of the curve represent the number of",
+        "pairs within a range. Note that exclusions do not affect the",
+        "normalization: even if [TT]-excl[tt] is set, or [TT]-ref[tt] and",
+        "[TT]-sel[tt] contain the same selection, the normalization factor",
+        "is still N*M, not N*(M-excluded).[PAR]",
+        "For [TT]-surf[tt], the selection provided to [TT]-ref[tt] must",
+        "select atoms, i.e., centers of mass are not supported. Further,",
+        "[TT]-nonorm[tt] is implied, as the bins have irregular shapes and",
+        "the volume of a bin is not easily computable.[PAR]",
+        "Option [TT]-cn[tt] produces the cumulative number RDF,",
+        "i.e. the average number of particles within a distance r.[PAR]"
+    };
+
+    options->setDescription(desc);
+
+    options->addOption(FileNameOption("o").filetype(eftPlot).outputFile().required()
+                           .store(&fnRdf_).defaultBasename("rdf")
+                           .description("Computed RDFs"));
+    options->addOption(FileNameOption("cn").filetype(eftPlot).outputFile()
+                           .store(&fnCumulative_).defaultBasename("rdf_cn")
+                           .description("Cumulative RDFs"));
+
+    options->addOption(DoubleOption("bin").store(&binwidth_)
+                           .description("Bin width (nm)"));
+    options->addOption(BooleanOption("norm").store(&bNormalize_)
+                           .description("Normalize for bin volume and density"));
+    options->addOption(BooleanOption("xy").store(&bXY_)
+                           .description("Use only the x and y components of the distance"));
+    options->addOption(BooleanOption("excl").store(&bExclusions_)
+                           .description("Use exclusions from topology"));
+    options->addOption(DoubleOption("cut").store(&cutoff_)
+                           .description("Shortest distance (nm) to be considered"));
+    options->addOption(DoubleOption("rmax").store(&rmax_)
+                           .description("Largest distance (nm) to calculate"));
+
+    const char *const cSurfaceEnum[] = { "no", "mol", "res" };
+    options->addOption(StringOption("surf").enumValue(cSurfaceEnum)
+                           .defaultEnumIndex(0).store(&surface_)
+                           .description("RDF with respect to the surface of the reference"));
+
+    options->addOption(SelectionOption("ref").store(&refSel_).required()
+                           .description("Reference selection for RDF computation"));
+    options->addOption(SelectionOption("sel").storeVector(&sel_)
+                           .required().multiValue()
+                           .description("Selections to compute RDFs for from the reference"));
+}
+
+void
+Rdf::optionsFinished(Options *options, TrajectoryAnalysisSettings *settings)
+{
+    if (surface_ != "no")
+    {
+        settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
+
+        if (options->isSet("norm") && bNormalize_)
+        {
+            GMX_THROW(InconsistentInputError("-surf cannot be combined with -norm"));
+        }
+        bNormalize_ = false;
+        if (bExclusions_)
+        {
+            GMX_THROW(InconsistentInputError("-surf cannot be combined with -excl"));
+        }
+    }
+    if (bExclusions_)
+    {
+        settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
+    }
+    if (cutoff_ < 0.0)
+    {
+        cutoff_ = 0.0;
+    }
+}
+
+void
+Rdf::initAnalysis(const TrajectoryAnalysisSettings &settings,
+                  const TopologyInformation        &top)
+{
+    pairDist_.setDataSetCount(sel_.size());
+    for (size_t i = 0; i < sel_.size(); ++i)
+    {
+        pairDist_.setColumnCount(i, 1);
+    }
+    plotSettings_ = settings.plotSettings();
+    nb_.setXYMode(bXY_);
+
+    normFactors_.setColumnCount(0, sel_.size() + 1);
+
+    const bool bSurface = (surface_ != "no");
+    if (bSurface)
+    {
+        if (!refSel_.hasOnlyAtoms())
+        {
+            GMX_THROW(InconsistentInputError("-surf only works with -refsel that consists of atoms"));
+        }
+        const e_index_t type = (surface_ == "mol" ? INDEX_MOL : INDEX_RES);
+        surfaceGroupCount_ = refSel_.initOriginalIdsToGroup(top.topology(), type);
+    }
+
+    if (bExclusions_)
+    {
+        if (!refSel_.hasOnlyAtoms())
+        {
+            GMX_THROW(InconsistentInputError("-excl only works with selections that consist of atoms"));
+        }
+        for (size_t i = 0; i < sel_.size(); ++i)
+        {
+            if (!sel_[i].hasOnlyAtoms())
+            {
+                GMX_THROW(InconsistentInputError("-excl only works with selections that consist of atoms"));
+            }
+        }
+        const t_topology *topology = top.topology();
+        if (topology->excls.nr == 0)
+        {
+            GMX_THROW(InconsistentInputError("-excl is set, but the file provided to -s does not define exclusions"));
+        }
+        nb_.setTopologyExclusions(&topology->excls);
+    }
+}
+
+void
+Rdf::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
+                         const t_trxframe                 &fr)
+{
+    // If -rmax is not provided, determine one from the box for the first frame.
+    if (rmax_ <= 0.0)
+    {
+        matrix box;
+        copy_mat(fr.box, box);
+        if (settings.hasPBC())
+        {
+            if (bXY_)
+            {
+                box[ZZ][ZZ] = 2*std::max(box[XX][XX], box[YY][YY]);
+            }
+            rmax_ = std::sqrt(0.99*0.99*max_cutoff2(bXY_ ? epbcXY : epbcXYZ, box));
+        }
+        else
+        {
+            if (bXY_)
+            {
+                clear_rvec(box[ZZ]);
+            }
+            rmax_ = 3*std::max(box[XX][XX], std::max(box[YY][YY], box[ZZ][ZZ]));
+        }
+    }
+    cut2_  = sqr(cutoff_);
+    rmax2_ = sqr(rmax_);
+    nb_.setCutoff(rmax_);
+    // We use the double amount of bins, so we can correctly
+    // write the rdf and rdf_cn output at i*binwidth values.
+    pairCounts_->init(histogramFromRange(0.0, rmax_).binWidth(binwidth_ / 2.0));
+}
+
+/*! \brief
+ * Temporary memory for use within a single-frame calculation.
+ */
+class RdfModuleData : public TrajectoryAnalysisModuleData
+{
+    public:
+        /*! \brief
+         * Reserves memory for the frame-local data.
+         *
+         * `surfaceGroupCount` will be zero if -surf is not specified.
+         */
+        RdfModuleData(TrajectoryAnalysisModule          *module,
+                      const AnalysisDataParallelOptions &opt,
+                      const SelectionCollection         &selections,
+                      int                                surfaceGroupCount)
+            : TrajectoryAnalysisModuleData(module, opt, selections)
+        {
+            surfaceDist2_.resize(surfaceGroupCount);
+        }
+
+        virtual void finish() { finishDataHandles(); }
+
+        /*! \brief
+         * Minimum distance to each surface group.
+         *
+         * One entry for each group (residue/molecule, per -surf) in the
+         * reference selection.
+         * This is needed to support neighborhood searching, which may not
+         * return the reference positions in order: for each position, we need
+         * to search through all the reference positions and update this array
+         * to find the minimum distance to each surface group, and then compute
+         * the RDF from these numbers.
+         */
+        std::vector<real> surfaceDist2_;
+};
+
+TrajectoryAnalysisModuleDataPointer Rdf::startFrames(
+        const AnalysisDataParallelOptions &opt,
+        const SelectionCollection         &selections)
+{
+    return TrajectoryAnalysisModuleDataPointer(
+            new RdfModuleData(this, opt, selections, surfaceGroupCount_));
+}
+
+void
+Rdf::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
+                  TrajectoryAnalysisModuleData *pdata)
+{
+    AnalysisDataHandle   dh        = pdata->dataHandle(pairDist_);
+    AnalysisDataHandle   nh        = pdata->dataHandle(normFactors_);
+    const Selection     &refSel    = pdata->parallelSelection(refSel_);
+    const SelectionList &sel       = pdata->parallelSelections(sel_);
+    RdfModuleData       &frameData = *static_cast<RdfModuleData *>(pdata);
+    const bool           bSurface  = !frameData.surfaceDist2_.empty();
+
+    matrix               boxForVolume;
+    copy_mat(fr.box, boxForVolume);
+    if (bXY_)
+    {
+        // Set z-size to 1 so we get the surface are iso the volume
+        clear_rvec(boxForVolume[ZZ]);
+        boxForVolume[ZZ][ZZ] = 1;
+    }
+    const real inverseVolume = 1.0 / det(boxForVolume);
+
+    nh.startFrame(frnr, fr.time);
+    // Compute the normalization factor for the number of reference positions.
+    if (bSurface)
+    {
+        if (refSel.isDynamic())
+        {
+            // Count the number of distinct groups.
+            // This assumes that each group is continuous, which is currently
+            // the case.
+            int count  = 0;
+            int prevId = -1;
+            for (int i = 0; i < refSel.posCount(); ++i)
+            {
+                const int id = refSel.position(i).mappedId();
+                if (id != prevId)
+                {
+                    ++count;
+                    prevId = id;
+                }
+            }
+            nh.setPoint(0, count);
+        }
+        else
+        {
+            nh.setPoint(0, surfaceGroupCount_);
+        }
+    }
+    else
+    {
+        nh.setPoint(0, refSel.posCount());
+    }
+
+    dh.startFrame(frnr, fr.time);
+    AnalysisNeighborhoodSearch    nbsearch = nb_.initSearch(pbc, refSel);
+    for (size_t g = 0; g < sel.size(); ++g)
+    {
+        dh.selectDataSet(g);
+
+        if (bSurface)
+        {
+            // Special loop for surface calculation, where a separate neighbor
+            // search is done for each position in the selection, and the
+            // nearest position from each surface group is tracked.
+            std::vector<real> &surfaceDist2 = frameData.surfaceDist2_;
+            for (int i = 0; i < sel[g].posCount(); ++i)
+            {
+                std::fill(surfaceDist2.begin(), surfaceDist2.end(),
+                          std::numeric_limits<real>::max());
+                AnalysisNeighborhoodPairSearch pairSearch =
+                    nbsearch.startPairSearch(sel[g].position(i));
+                AnalysisNeighborhoodPair       pair;
+                while (pairSearch.findNextPair(&pair))
+                {
+                    const real r2    = pair.distance2();
+                    const int  refId = refSel.position(pair.refIndex()).mappedId();
+                    if (r2 < surfaceDist2[refId])
+                    {
+                        surfaceDist2[refId] = r2;
+                    }
+                }
+                // Accumulate the RDF from the distances to the surface.
+                for (size_t i = 0; i < surfaceDist2.size(); ++i)
+                {
+                    const real r2 = surfaceDist2[i];
+                    // Here, we need to check for rmax, since the value might
+                    // be above the cutoff if no points were close to some
+                    // surface positions.
+                    if (r2 > cut2_ && r2 <= rmax2_)
+                    {
+                        dh.setPoint(0, std::sqrt(r2));
+                        dh.finishPointSet();
+                    }
+                }
+            }
+        }
+        else
+        {
+            // Standard neighborhood search over all pairs within the cutoff
+            // for the -surf no case.
+            AnalysisNeighborhoodPairSearch pairSearch = nbsearch.startPairSearch(sel[g]);
+            AnalysisNeighborhoodPair       pair;
+            while (pairSearch.findNextPair(&pair))
+            {
+                const real r2 = pair.distance2();
+                if (r2 > cut2_)
+                {
+                    // TODO: Consider whether the histogramming could be done with
+                    // less overhead (after first measuring the overhead).
+                    dh.setPoint(0, std::sqrt(r2));
+                    dh.finishPointSet();
+                }
+            }
+        }
+        // Normalization factor for the number density (only used without
+        // -surf, but does not hurt to populate otherwise).
+        nh.setPoint(g + 1, sel[g].posCount() * inverseVolume);
+    }
+    dh.finishFrame();
+    nh.finishFrame();
+}
+
+void
+Rdf::finishAnalysis(int /*nframes*/)
+{
+    // Normalize the averager with the number of reference positions,
+    // from where the normalization propagates to all the output.
+    const real refPosCount = normAve_->average(0, 0);
+    pairCounts_->averager().scaleAll(1.0 / refPosCount);
+    pairCounts_->averager().done();
+
+    // TODO: Consider how these could be exposed to the testing framework
+    // through the dataset registration mechanism.
+    AverageHistogramPointer finalRdf
+        = pairCounts_->averager().resampleDoubleBinWidth(true);
+
+    if (bNormalize_)
+    {
+        // Normalize by the volume of the bins (volume of sphere segments or
+        // length of circle segments).
+        std::vector<real> invBinVolume;
+        const int         nbin = finalRdf->settings().binCount();
+        invBinVolume.resize(nbin);
+        real              prevSphereVolume = 0.0;
+        for (int i = 0; i < nbin; ++i)
+        {
+            const real r = (i + 0.5)*binwidth_;
+            real       sphereVolume;
+            if (bXY_)
+            {
+                sphereVolume = M_PI*r*r;
+            }
+            else
+            {
+                sphereVolume = (4.0/3.0)*M_PI*r*r*r;
+            }
+            const real binVolume = sphereVolume - prevSphereVolume;
+            invBinVolume[i]  = 1.0 / binVolume;
+            prevSphereVolume = sphereVolume;
+        }
+        finalRdf->scaleAllByVector(&invBinVolume[0]);
+
+        // Normalize by particle density.
+        for (size_t g = 0; g < sel_.size(); ++g)
+        {
+            finalRdf->scaleSingle(g, 1.0 / normAve_->average(0, g + 1));
+        }
+    }
+    else
+    {
+        // With -nonorm, just scale with bin width to make the integral of the
+        // curve (instead of raw bin sum) represent the pair count.
+        finalRdf->scaleAll(1.0 / binwidth_);
+    }
+    finalRdf->done();
+
+    // TODO: Consider if some of this should be done in writeOutput().
+    {
+        AnalysisDataPlotModulePointer plotm(
+                new AnalysisDataPlotModule(plotSettings_));
+        plotm->setFileName(fnRdf_);
+        plotm->setTitle("Radial distribution");
+        plotm->setSubtitle(formatString("reference %s", refSel_.name()));
+        plotm->setXLabel("r (nm)");
+        plotm->setYLabel("g(r)");
+        for (size_t i = 0; i < sel_.size(); ++i)
+        {
+            plotm->appendLegend(sel_[i].name());
+        }
+        finalRdf->addModule(plotm);
+    }
+
+    if (!fnCumulative_.empty())
+    {
+        AverageHistogramPointer cumulativeRdf
+            = pairCounts_->averager().resampleDoubleBinWidth(false);
+        cumulativeRdf->makeCumulative();
+        cumulativeRdf->done();
+
+        AnalysisDataPlotModulePointer plotm(
+                new AnalysisDataPlotModule(plotSettings_));
+        plotm->setFileName(fnCumulative_);
+        plotm->setTitle("Cumulative Number RDF");
+        plotm->setSubtitle(formatString("reference %s", refSel_.name()));
+        plotm->setXLabel("r (nm)");
+        plotm->setYLabel("number");
+        for (size_t i = 0; i < sel_.size(); ++i)
+        {
+            plotm->appendLegend(sel_[i].name());
+        }
+        cumulativeRdf->addModule(plotm);
+    }
+}
+
+void
+Rdf::writeOutput()
+{
+}
+
+//! \}
+
+}       // namespace
+
+const char RdfInfo::name[]             = "rdf";
+const char RdfInfo::shortDescription[] =
+    "Calculate radial distribution functions";
+
+TrajectoryAnalysisModulePointer RdfInfo::create()
+{
+    return TrajectoryAnalysisModulePointer(new Rdf);
+}
+
+} // namespace analysismodules
+
+} // namespace gmx
diff --git a/src/gromacs/trajectoryanalysis/modules/rdf.h b/src/gromacs/trajectoryanalysis/modules/rdf.h
new file mode 100644 (file)
index 0000000..269d2a8
--- /dev/null
@@ -0,0 +1,65 @@
+/*
+ * 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
+ * Declares trajectory analysis module for RDF calculations.
+ *
+ * \author Teemu Murtola <teemu.murtola@gmail.com>
+ * \ingroup module_trajectoryanalysis
+ */
+#ifndef GMX_TRAJECTORYANALYSIS_MODULES_RDF_H
+#define GMX_TRAJECTORYANALYSIS_MODULES_RDF_H
+
+#include "gromacs/trajectoryanalysis/analysismodule.h"
+
+namespace gmx
+{
+
+namespace analysismodules
+{
+
+class RdfInfo
+{
+    public:
+        static const char name[];
+        static const char shortDescription[];
+        static TrajectoryAnalysisModulePointer create();
+};
+
+} // namespace analysismodules
+
+} // namespace gmx
+
+#endif
index 4fafb502e41a0df02646df1bf2b096e8de169452..27bfbd17c4a8cbbf0ac61bb30a76e9b42dfc5dd2 100644 (file)
@@ -38,6 +38,7 @@ gmx_add_unit_test(TrajectoryAnalysisUnitTests trajectoryanalysis-test
                   distance.cpp
                   freevolume.cpp
                   pairdist.cpp
+                  rdf.cpp
                   sasa.cpp
                   select.cpp
                   surfacearea.cpp
diff --git a/src/gromacs/trajectoryanalysis/tests/rdf.cpp b/src/gromacs/trajectoryanalysis/tests/rdf.cpp
new file mode 100644 (file)
index 0000000..db97f2d
--- /dev/null
@@ -0,0 +1,120 @@
+/*
+ * 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 functionality of the "rdf" trajectory analysis module.
+ *
+ * These tests are essentially regression tests for the actual RDF calculation
+ * from very small configurations.  Exclusions are not tested (since the input
+ * does not contain any), nor is the effect of -cut of -rmax.
+ * At the moment, they do not test the final normalization, but only the pair
+ * counts calculated for each frame.  Tests for the final normalization should
+ * be added once related TODOs in the implementation/framework have been
+ * resolved (currently, the test framework does not see the final output data).
+ *
+ * \author Teemu Murtola <teemu.murtola@gmail.com>
+ * \ingroup module_trajectoryanalysis
+ */
+#include "gmxpre.h"
+
+#include "gromacs/trajectoryanalysis/modules/rdf.h"
+
+#include <gtest/gtest.h>
+
+#include "testutils/cmdlinetest.h"
+#include "testutils/testasserts.h"
+
+#include "moduletest.h"
+
+namespace
+{
+
+using gmx::test::CommandLine;
+
+/********************************************************************
+ * Tests for gmx::analysismodules::Rdf.
+ */
+
+//! Test fixture for the `rdf` analysis module.
+typedef gmx::test::TrajectoryAnalysisModuleTestFixture<gmx::analysismodules::RdfInfo>
+    RdfModuleTest;
+
+TEST_F(RdfModuleTest, BasicTest)
+{
+    const char *const cmdline[] = {
+        "rdf",
+        "-bin", "0.05",
+        "-ref", "name OW",
+        "-sel", "name OW", "not name OW"
+    };
+    setTopology("spc216.gro");
+    setOutputFileNoTest("-o", "xvg");
+    excludeDataset("pairdist");
+    runTest(CommandLine(cmdline));
+}
+
+TEST_F(RdfModuleTest, CalculatesSurf)
+{
+    const char *const cmdline[] = {
+        "rdf",
+        "-bin", "0.05", "-surf", "res",
+        "-ref", "within 0.5 of (resnr 1 and name OW)",
+        "-sel", "name OW", "not name OW"
+    };
+    setTopology("spc216.gro");
+    setOutputFileNoTest("-o", "xvg");
+    excludeDataset("pairdist");
+    runTest(CommandLine(cmdline));
+}
+
+TEST_F(RdfModuleTest, CalculatesXY)
+{
+    const char *const cmdline[] = {
+        "rdf",
+        "-bin", "0.05", "-xy",
+        "-ref", "name OW",
+        "-sel", "name OW", "not name OW"
+    };
+    setTopology("spc216.gro");
+    setOutputFileNoTest("-o", "xvg");
+    excludeDataset("pairdist");
+    // TODO: Consider if it is possible to get a more reproducible result
+    // and/or a stricter tolerance (e.g., by checking that the sum of
+    // neighboring values still stays constant).
+    setDatasetTolerance("paircount", gmx::test::absoluteTolerance(2.5));
+    runTest(CommandLine(cmdline));
+}
+
+} // namespace
diff --git a/src/gromacs/trajectoryanalysis/tests/refdata/RdfModuleTest_BasicTest.xml b/src/gromacs/trajectoryanalysis/tests/refdata/RdfModuleTest_BasicTest.xml
new file mode 100644 (file)
index 0000000..8fff980
--- /dev/null
@@ -0,0 +1,259 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <String Name="CommandLine">rdf -bin 0.05 -ref 'name OW' -sel 'name OW' 'not name OW'</String>
+  <OutputData Name="Data">
+    <AnalysisData Name="norm">
+      <DataFrame Name="Frame0">
+        <Real Name="X">0</Real>
+        <DataValues>
+          <Int Name="Count">3</Int>
+          <DataValue>
+            <Real Name="Value">216</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">33.455901630929986</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">66.911803261859973</Real>
+          </DataValue>
+        </DataValues>
+      </DataFrame>
+    </AnalysisData>
+    <AnalysisData Name="paircount">
+      <DataFrame Name="Frame0">
+        <Real Name="X">0</Real>
+        <DataValues>
+          <Int Name="Count">37</Int>
+          <Int Name="DataSet">0</Int>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">274</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">360</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">226</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">234</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">270</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">332</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">420</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">456</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">548</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">588</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">546</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">632</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">660</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">696</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">822</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">922</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1060</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1084</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1276</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1260</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1260</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1416</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1468</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1560</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1668</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1774</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1578</Real>
+          </DataValue>
+        </DataValues>
+        <DataValues>
+          <Int Name="Count">37</Int>
+          <Int Name="DataSet">1</Int>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">215</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">217</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">114</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">163</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">87</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">52</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">103</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">266</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">618</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">751</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">703</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">722</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">772</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">821</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">946</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1065</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1229</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1281</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1402</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1518</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1640</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1844</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2058</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2137</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2412</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2451</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2650</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2886</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2928</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">3101</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">3374</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">3623</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">3288</Real>
+          </DataValue>
+        </DataValues>
+      </DataFrame>
+    </AnalysisData>
+  </OutputData>
+</ReferenceData>
diff --git a/src/gromacs/trajectoryanalysis/tests/refdata/RdfModuleTest_CalculatesSurf.xml b/src/gromacs/trajectoryanalysis/tests/refdata/RdfModuleTest_CalculatesSurf.xml
new file mode 100644 (file)
index 0000000..1c8be52
--- /dev/null
@@ -0,0 +1,259 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <String Name="CommandLine">rdf -bin 0.05 -surf res -ref 'within 0.5 of (resnr 1 and name OW)' -sel 'name OW' 'not name OW'</String>
+  <OutputData Name="Data">
+    <AnalysisData Name="norm">
+      <DataFrame Name="Frame0">
+        <Real Name="X">0</Real>
+        <DataValues>
+          <Int Name="Count">3</Int>
+          <DataValue>
+            <Real Name="Value">22</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">33.455901630929986</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">66.911803261859973</Real>
+          </DataValue>
+        </DataValues>
+      </DataFrame>
+    </AnalysisData>
+    <AnalysisData Name="paircount">
+      <DataFrame Name="Frame0">
+        <Real Name="X">0</Real>
+        <DataValues>
+          <Int Name="Count">37</Int>
+          <Int Name="DataSet">0</Int>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">7</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">15</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">8</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">3</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">22</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">28</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">28</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">42</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">38</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">41</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">52</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">48</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">57</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">77</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">58</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">80</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">83</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">82</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">90</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">98</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">128</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">127</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">136</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">130</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">155</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">151</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">189</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">153</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">183</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">201</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">153</Real>
+          </DataValue>
+        </DataValues>
+        <DataValues>
+          <Int Name="Count">37</Int>
+          <Int Name="DataSet">1</Int>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">3</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">10</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">22</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">21</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">33</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">51</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">43</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">61</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">74</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">61</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">80</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">86</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">116</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">136</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">138</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">123</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">137</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">177</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">168</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">179</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">205</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">241</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">239</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">303</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">261</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">303</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">319</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">321</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">354</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">362</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">374</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">352</Real>
+          </DataValue>
+        </DataValues>
+      </DataFrame>
+    </AnalysisData>
+  </OutputData>
+</ReferenceData>
diff --git a/src/gromacs/trajectoryanalysis/tests/refdata/RdfModuleTest_CalculatesXY.xml b/src/gromacs/trajectoryanalysis/tests/refdata/RdfModuleTest_CalculatesXY.xml
new file mode 100644 (file)
index 0000000..6fb0ae3
--- /dev/null
@@ -0,0 +1,259 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <String Name="CommandLine">rdf -bin 0.05 -xy -ref 'name OW' -sel 'name OW' 'not name OW'</String>
+  <OutputData Name="Data">
+    <AnalysisData Name="norm">
+      <DataFrame Name="Frame0">
+        <Real Name="X">0</Real>
+        <DataValues>
+          <Int Name="Count">3</Int>
+          <DataValue>
+            <Real Name="Value">216</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">62.296896190889484</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">124.59379238177897</Real>
+          </DataValue>
+        </DataValues>
+      </DataFrame>
+    </AnalysisData>
+    <AnalysisData Name="paircount">
+      <DataFrame Name="Frame0">
+        <Real Name="X">0</Real>
+        <DataValues>
+          <Int Name="Count">37</Int>
+          <Int Name="DataSet">0</Int>
+          <DataValue>
+            <Real Name="Value">28</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">66</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">90</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">150</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">206</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">222</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">276</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">366</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">420</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">522</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">600</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">632</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">644</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">766</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">820</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">818</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">828</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">932</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">950</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1042</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1088</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1030</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1184</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1276</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1346</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1388</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1380</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1466</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1572</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1474</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1578</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1630</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1782</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1754</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1772</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1854</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1686</Real>
+          </DataValue>
+        </DataValues>
+        <DataValues>
+          <Int Name="Count">37</Int>
+          <Int Name="DataSet">1</Int>
+          <DataValue>
+            <Real Name="Value">64</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">184</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">288</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">567</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">408</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">539</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">630</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">743</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">821</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">897</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1094</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1236</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1394</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1504</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1568</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1691</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1770</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1827</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1907</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2001</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2170</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2194</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2367</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2521</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2546</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2731</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2909</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">2962</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">3044</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">3038</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">3139</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">3296</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">3507</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">3404</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">3612</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">3945</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">3388</Real>
+          </DataValue>
+        </DataValues>
+      </DataFrame>
+    </AnalysisData>
+  </OutputData>
+</ReferenceData>
index 55515b411568350738622349ff21b4da0d0323d4..ca387efd4ceed24b4febbf5111dccd356f39f152 100644 (file)
@@ -68,10 +68,10 @@ namespace
 using gmx::test::CommandLine;
 
 /********************************************************************
- * Tests for gmx::analysismodules::Sas.
+ * Tests for gmx::analysismodules::Sasa.
  */
 
-//! Test fixture for the select analysis module.
+//! Test fixture for the `sasa` analysis module.
 typedef gmx::test::TrajectoryAnalysisModuleTestFixture<gmx::analysismodules::SasaInfo>
     SasaModuleTest;
 
diff --git a/src/gromacs/trajectoryanalysis/tests/spc216.gro b/src/gromacs/trajectoryanalysis/tests/spc216.gro
new file mode 100644 (file)
index 0000000..49449af
--- /dev/null
@@ -0,0 +1,651 @@
+216H2O,WATJP01,SPC216,SPC-MODEL,300K,BOX(M)=1.86206NM,WFVG,MAR. 1984
+  648
+    1SOL     OW    1    .230    .628    .113
+    1SOL    HW1    2    .137    .626    .150
+    1SOL    HW2    3    .231    .589    .021
+    2SOL     OW    4    .225    .275   -.866
+    2SOL    HW1    5    .260    .258   -.774
+    2SOL    HW2    6    .137    .230   -.878
+    3SOL     OW    7    .019    .368    .647
+    3SOL    HW1    8   -.063    .411    .686
+    3SOL    HW2    9   -.009    .295    .584
+    4SOL     OW   10    .569   -.587   -.697
+    4SOL    HW1   11    .476   -.594   -.734
+    4SOL    HW2   12    .580   -.498   -.653
+    5SOL     OW   13   -.307   -.351    .703
+    5SOL    HW1   14   -.364   -.367    .784
+    5SOL    HW2   15   -.366   -.341    .623
+    6SOL     OW   16   -.119    .618    .856
+    6SOL    HW1   17   -.086    .712    .856
+    6SOL    HW2   18   -.068    .564    .922
+    7SOL     OW   19   -.727    .703    .717
+    7SOL    HW1   20   -.670    .781    .692
+    7SOL    HW2   21   -.787    .729    .793
+    8SOL     OW   22   -.107    .607    .231
+    8SOL    HW1   23   -.119    .594    .132
+    8SOL    HW2   24   -.137    .526    .280
+    9SOL     OW   25    .768   -.718   -.839
+    9SOL    HW1   26    .690   -.701   -.779
+    9SOL    HW2   27    .802   -.631   -.875
+   10SOL     OW   28    .850    .798   -.039
+   10SOL    HW1   29    .846    .874    .026
+   10SOL    HW2   30    .872    .834   -.130
+   11SOL     OW   31    .685   -.850    .665
+   11SOL    HW1   32    .754   -.866    .735
+   11SOL    HW2   33    .612   -.793    .703
+   12SOL     OW   34    .686   -.701   -.059
+   12SOL    HW1   35    .746   -.622   -.045
+   12SOL    HW2   36    .600   -.670   -.100
+   13SOL     OW   37    .335   -.427   -.801
+   13SOL    HW1   38    .257   -.458   -.854
+   13SOL    HW2   39    .393   -.369   -.858
+   14SOL     OW   40   -.402   -.357   -.523
+   14SOL    HW1   41   -.378   -.263   -.497
+   14SOL    HW2   42   -.418   -.411   -.441
+   15SOL     OW   43    .438    .392   -.363
+   15SOL    HW1   44    .520    .336   -.354
+   15SOL    HW2   45    .357    .334   -.359
+   16SOL     OW   46   -.259    .447    .737
+   16SOL    HW1   47   -.333    .493    .687
+   16SOL    HW2   48   -.208    .515    .790
+   17SOL     OW   49    .231   -.149    .483
+   17SOL    HW1   50    .265   -.072    .537
+   17SOL    HW2   51    .275   -.149    .393
+   18SOL     OW   52   -.735   -.521   -.172
+   18SOL    HW1   53   -.688   -.521   -.084
+   18SOL    HW2   54   -.783   -.608   -.183
+   19SOL     OW   55    .230   -.428    .538
+   19SOL    HW1   56    .204   -.332    .538
+   19SOL    HW2   57    .159   -.482    .583
+   20SOL     OW   58    .240   -.771    .886
+   20SOL    HW1   59    .254   -.855    .938
+   20SOL    HW2   60    .185   -.707    .941
+   21SOL     OW   61    .620   -.076   -.423
+   21SOL    HW1   62    .528   -.093   -.388
+   21SOL    HW2   63    .648    .016   -.397
+   22SOL     OW   64    .606   -.898    .123
+   22SOL    HW1   65    .613   -.814    .069
+   22SOL    HW2   66    .652   -.885    .211
+   23SOL     OW   67   -.268    .114   -.382
+   23SOL    HW1   68   -.286    .181   -.454
+   23SOL    HW2   69   -.271    .160   -.293
+   24SOL     OW   70    .122    .643    .563
+   24SOL    HW1   71    .077    .555    .580
+   24SOL    HW2   72    .121    .697    .647
+   25SOL     OW   73   -.020   -.095    .359
+   25SOL    HW1   74    .034   -.124    .439
+   25SOL    HW2   75    .010   -.005    .330
+   26SOL     OW   76    .027   -.266    .117
+   26SOL    HW1   77    .008   -.362    .138
+   26SOL    HW2   78   -.006   -.208    .192
+   27SOL     OW   79   -.173    .922    .612
+   27SOL    HW1   80   -.078    .893    .620
+   27SOL    HW2   81   -.181    .987    .537
+   28SOL     OW   82   -.221   -.754    .432
+   28SOL    HW1   83   -.135   -.752    .380
+   28SOL    HW2   84   -.207   -.707    .520
+   29SOL     OW   85    .113    .737   -.265
+   29SOL    HW1   86    .201    .724   -.220
+   29SOL    HW2   87    .100    .834   -.287
+   30SOL     OW   88    .613   -.497    .726
+   30SOL    HW1   89    .564   -.584    .735
+   30SOL    HW2   90    .590   -.454    .639
+   31SOL     OW   91   -.569   -.634   -.439
+   31SOL    HW1   92   -.532   -.707   -.497
+   31SOL    HW2   93   -.517   -.629   -.354
+   32SOL     OW   94    .809    .004    .502
+   32SOL    HW1   95    .849    .095    .493
+   32SOL    HW2   96    .709    .012    .508
+   33SOL     OW   97    .197   -.886   -.598
+   33SOL    HW1   98    .286   -.931   -.612
+   33SOL    HW2   99    .124   -.951   -.617
+   34SOL     OW  100   -.337   -.863    .190
+   34SOL    HW1  101   -.400   -.939    .203
+   34SOL    HW2  102   -.289   -.845    .276
+   35SOL     OW  103   -.675   -.070   -.246
+   35SOL    HW1  104   -.651   -.010   -.322
+   35SOL    HW2  105   -.668   -.165   -.276
+   36SOL     OW  106    .317    .251   -.061
+   36SOL    HW1  107    .388    .322   -.055
+   36SOL    HW2  108    .229    .290   -.033
+   37SOL     OW  109   -.396   -.445   -.909
+   37SOL    HW1  110   -.455   -.439   -.829
+   37SOL    HW2  111   -.411   -.533   -.955
+   38SOL     OW  112   -.195   -.148    .572
+   38SOL    HW1  113   -.236   -.171    .484
+   38SOL    HW2  114   -.213   -.222    .637
+   39SOL     OW  115    .598    .729    .270
+   39SOL    HW1  116    .622    .798    .202
+   39SOL    HW2  117    .520    .762    .324
+   40SOL     OW  118   -.581    .345   -.918
+   40SOL    HW1  119   -.667    .295   -.931
+   40SOL    HW2  120   -.519    .291   -.862
+   41SOL     OW  121   -.286   -.200    .307
+   41SOL    HW1  122   -.197   -.154    .310
+   41SOL    HW2  123   -.307   -.224    .212
+   42SOL     OW  124    .807    .605   -.397
+   42SOL    HW1  125    .760    .602   -.308
+   42SOL    HW2  126    .756    .550   -.463
+   43SOL     OW  127   -.468    .469   -.188
+   43SOL    HW1  128   -.488    .512   -.100
+   43SOL    HW2  129   -.390    .407   -.179
+   44SOL     OW  130   -.889    .890   -.290
+   44SOL    HW1  131   -.843    .806   -.319
+   44SOL    HW2  132   -.945    .924   -.365
+   45SOL     OW  133   -.871    .410   -.620
+   45SOL    HW1  134   -.948    .444   -.566
+   45SOL    HW2  135   -.905    .359   -.699
+   46SOL     OW  136   -.821    .701    .429
+   46SOL    HW1  137   -.795    .697    .525
+   46SOL    HW2  138   -.906    .650    .415
+   47SOL     OW  139    .076    .811    .789
+   47SOL    HW1  140    .175    .799    .798
+   47SOL    HW2  141    .052    .906    .810
+   48SOL     OW  142    .130   -.041   -.291
+   48SOL    HW1  143    .120   -.056   -.192
+   48SOL    HW2  144    .044   -.005   -.327
+   49SOL     OW  145    .865    .348    .195
+   49SOL    HW1  146    .924    .411    .146
+   49SOL    HW2  147    .884    .254    .166
+   50SOL     OW  148   -.143    .585   -.031
+   50SOL    HW1  149   -.169    .674   -.067
+   50SOL    HW2  150   -.145    .517   -.104
+   51SOL     OW  151   -.500   -.718    .545
+   51SOL    HW1  152   -.417   -.747    .497
+   51SOL    HW2  153   -.549   -.651    .489
+   52SOL     OW  154    .550    .196    .885
+   52SOL    HW1  155    .545    .191    .985
+   52SOL    HW2  156    .552    .292    .856
+   53SOL     OW  157   -.854   -.406    .477
+   53SOL    HW1  158   -.900   -.334    .425
+   53SOL    HW2  159   -.858   -.386    .575
+   54SOL     OW  160    .351   -.061    .853
+   54SOL    HW1  161    .401   -.147    .859
+   54SOL    HW2  162    .416    .016    .850
+   55SOL     OW  163   -.067   -.796    .873
+   55SOL    HW1  164   -.129   -.811    .797
+   55SOL    HW2  165   -.119   -.785    .958
+   56SOL     OW  166   -.635   -.312   -.356
+   56SOL    HW1  167   -.629   -.389   -.292
+   56SOL    HW2  168   -.687   -.338   -.436
+   57SOL     OW  169    .321   -.919    .242
+   57SOL    HW1  170    .403   -.880    .200
+   57SOL    HW2  171    .294  -1.001    .193
+   58SOL     OW  172   -.404    .735    .728
+   58SOL    HW1  173   -.409    .670    .803
+   58SOL    HW2  174   -.324    .794    .741
+   59SOL     OW  175    .461   -.596   -.135
+   59SOL    HW1  176    .411   -.595   -.221
+   59SOL    HW2  177    .398   -.614   -.059
+   60SOL     OW  178   -.751   -.086    .237
+   60SOL    HW1  179   -.811   -.148    .287
+   60SOL    HW2  180   -.720   -.130    .152
+   61SOL     OW  181    .202    .285   -.364
+   61SOL    HW1  182    .122    .345   -.377
+   61SOL    HW2  183    .192    .236   -.278
+   62SOL     OW  184   -.230   -.485    .081
+   62SOL    HW1  185   -.262   -.391    .071
+   62SOL    HW2  186   -.306   -.548    .069
+   63SOL     OW  187    .464   -.119    .323
+   63SOL    HW1  188    .497   -.080    .409
+   63SOL    HW2  189    .540   -.126    .258
+   64SOL     OW  190   -.462    .107    .426
+   64SOL    HW1  191   -.486    .070    .336
+   64SOL    HW2  192   -.363    .123    .430
+   65SOL     OW  193    .249   -.077   -.621
+   65SOL    HW1  194    .306   -.142   -.571
+   65SOL    HW2  195    .233   -.110   -.714
+   66SOL     OW  196   -.922   -.164    .904
+   66SOL    HW1  197   -.842   -.221    .925
+   66SOL    HW2  198   -.971   -.204    .827
+   67SOL     OW  199    .382    .700    .480
+   67SOL    HW1  200    .427    .610    .477
+   67SOL    HW2  201    .288    .689    .513
+   68SOL     OW  202   -.315    .222   -.133
+   68SOL    HW1  203   -.320    .259   -.041
+   68SOL    HW2  204   -.387    .153   -.145
+   69SOL     OW  205    .614    .122    .117
+   69SOL    HW1  206    .712    .100    .124
+   69SOL    HW2  207    .583    .105    .024
+   70SOL     OW  208    .781    .264   -.113
+   70SOL    HW1  209    .848    .203   -.070
+   70SOL    HW2  210    .708    .283   -.048
+   71SOL     OW  211    .888   -.348   -.667
+   71SOL    HW1  212    .865   -.373   -.761
+   71SOL    HW2  213    .949   -.417   -.628
+   72SOL     OW  214   -.511    .590   -.429
+   72SOL    HW1  215   -.483    .547   -.344
+   72SOL    HW2  216   -.486    .686   -.428
+   73SOL     OW  217    .803   -.460    .924
+   73SOL    HW1  218    .893   -.446    .882
+   73SOL    HW2  219    .732   -.458    .853
+   74SOL     OW  220    .922    .503    .899
+   74SOL    HW1  221    .897    .494    .803
+   74SOL    HW2  222    .970    .421    .930
+   75SOL     OW  223    .539    .064    .512
+   75SOL    HW1  224    .458    .065    .570
+   75SOL    HW2  225    .542    .147    .457
+   76SOL     OW  226   -.428   -.674    .041
+   76SOL    HW1  227   -.396   -.750    .098
+   76SOL    HW2  228   -.520   -.647    .071
+   77SOL     OW  229    .297    .035    .171
+   77SOL    HW1  230    .346    .119    .150
+   77SOL    HW2  231    .359   -.030    .216
+   78SOL     OW  232   -.927    .236    .480
+   78SOL    HW1  233   -.975    .277    .402
+   78SOL    HW2  234   -.828    .234    .461
+   79SOL     OW  235   -.786    .683   -.398
+   79SOL    HW1  236   -.866    .622   -.395
+   79SOL    HW2  237   -.705    .630   -.422
+   80SOL     OW  238   -.635   -.292    .793
+   80SOL    HW1  239   -.614   -.218    .728
+   80SOL    HW2  240   -.567   -.292    .866
+   81SOL     OW  241    .459   -.710    .741
+   81SOL    HW1  242    .388   -.737    .806
+   81SOL    HW2  243    .433   -.738    .648
+   82SOL     OW  244   -.591   -.065    .591
+   82SOL    HW1  245   -.547   -.001    .527
+   82SOL    HW2  246   -.641   -.013    .661
+   83SOL     OW  247   -.830    .549    .016
+   83SOL    HW1  248   -.871    .631   -.023
+   83SOL    HW2  249   -.766    .575    .089
+   84SOL     OW  250    .078    .556   -.476
+   84SOL    HW1  251    .170    .555   -.517
+   84SOL    HW2  252    .072    .630   -.409
+   85SOL     OW  253    .561    .222   -.715
+   85SOL    HW1  254    .599    .138   -.678
+   85SOL    HW2  255    .473    .241   -.671
+   86SOL     OW  256    .866    .454    .642
+   86SOL    HW1  257    .834    .526    .580
+   86SOL    HW2  258    .890    .373    .589
+   87SOL     OW  259   -.845    .039    .753
+   87SOL    HW1  260   -.917    .044    .684
+   87SOL    HW2  261   -.869   -.030    .822
+   88SOL     OW  262   -.433   -.689    .867
+   88SOL    HW1  263   -.488   -.773    .860
+   88SOL    HW2  264   -.407   -.660    .775
+   89SOL     OW  265   -.396    .590   -.870
+   89SOL    HW1  266   -.426    .495   -.863
+   89SOL    HW2  267   -.323    .606   -.804
+   90SOL     OW  268   -.005    .833    .377
+   90SOL    HW1  269    .037    .769    .441
+   90SOL    HW2  270   -.043    .782    .299
+   91SOL     OW  271    .488   -.477    .174
+   91SOL    HW1  272    .401   -.492    .221
+   91SOL    HW2  273    .471   -.451    .079
+   92SOL     OW  274   -.198   -.582    .657
+   92SOL    HW1  275   -.099   -.574    .671
+   92SOL    HW2  276   -.243   -.498    .688
+   93SOL     OW  277   -.472    .575    .078
+   93SOL    HW1  278   -.526    .554    .159
+   93SOL    HW2  279   -.381    .534    .087
+   94SOL     OW  280    .527    .256    .328
+   94SOL    HW1  281    .554    .197    .253
+   94SOL    HW2  282    .527    .351    .297
+   95SOL     OW  283   -.108   -.639   -.274
+   95SOL    HW1  284   -.017   -.678   -.287
+   95SOL    HW2  285   -.100   -.543   -.250
+   96SOL     OW  286   -.798   -.515   -.522
+   96SOL    HW1  287   -.878   -.538   -.467
+   96SOL    HW2  288   -.715   -.541   -.473
+   97SOL     OW  289   -.270   -.233   -.237
+   97SOL    HW1  290   -.243   -.199   -.327
+   97SOL    HW2  291   -.191   -.271   -.191
+   98SOL     OW  292   -.751   -.667   -.762
+   98SOL    HW1  293   -.791   -.623   -.681
+   98SOL    HW2  294   -.792   -.630   -.845
+   99SOL     OW  295   -.224   -.763   -.783
+   99SOL    HW1  296   -.219   -.682   -.724
+   99SOL    HW2  297   -.310   -.761   -.834
+  100SOL     OW  298    .915    .089   -.460
+  100SOL    HW1  299    .940    .069   -.555
+  100SOL    HW2  300    .987    .145   -.418
+  101SOL     OW  301   -.882   -.746   -.143
+  101SOL    HW1  302   -.981   -.740   -.133
+  101SOL    HW2  303   -.859   -.826   -.199
+  102SOL     OW  304    .705   -.812    .368
+  102SOL    HW1  305    .691   -.805    .467
+  102SOL    HW2  306    .789   -.863    .350
+  103SOL     OW  307    .410    .813   -.611
+  103SOL    HW1  308    .496    .825   -.561
+  103SOL    HW2  309    .368    .726   -.584
+  104SOL     OW  310   -.588    .386   -.600
+  104SOL    HW1  311   -.567    .460   -.536
+  104SOL    HW2  312   -.677    .403   -.643
+  105SOL     OW  313    .064   -.298   -.531
+  105SOL    HW1  314    .018   -.216   -.565
+  105SOL    HW2  315    .162   -.279   -.522
+  106SOL     OW  316    .367   -.762    .501
+  106SOL    HW1  317    .360   -.679    .445
+  106SOL    HW2  318    .371   -.842    .441
+  107SOL     OW  319    .566    .537    .865
+  107SOL    HW1  320    .578    .603    .791
+  107SOL    HW2  321    .612    .571    .948
+  108SOL     OW  322   -.610   -.514    .388
+  108SOL    HW1  323   -.560   -.437    .428
+  108SOL    HW2  324   -.705   -.512    .420
+  109SOL     OW  325   -.590   -.417   -.720
+  109SOL    HW1  326   -.543   -.404   -.633
+  109SOL    HW2  327   -.656   -.491   -.711
+  110SOL     OW  328   -.280    .639    .472
+  110SOL    HW1  329   -.311    .700    .545
+  110SOL    HW2  330   -.230    .691    .403
+  111SOL     OW  331    .354   -.352   -.533
+  111SOL    HW1  332    .333   -.396   -.620
+  111SOL    HW2  333    .451   -.326   -.530
+  112SOL     OW  334    .402    .751   -.264
+  112SOL    HW1  335    .470    .806   -.311
+  112SOL    HW2  336    .442    .663   -.237
+  113SOL     OW  337   -.275    .779   -.192
+  113SOL    HW1  338   -.367    .817   -.197
+  113SOL    HW2  339   -.215    .826   -.257
+  114SOL     OW  340   -.849    .105   -.092
+  114SOL    HW1  341   -.843    .190   -.144
+  114SOL    HW2  342   -.817    .029   -.149
+  115SOL     OW  343    .504    .050   -.122
+  115SOL    HW1  344    .462   -.007   -.192
+  115SOL    HW2  345    .438    .119   -.090
+  116SOL     OW  346    .573    .870   -.833
+  116SOL    HW1  347    .617    .959   -.842
+  116SOL    HW2  348    .510    .870   -.756
+  117SOL     OW  349   -.502    .862   -.817
+  117SOL    HW1  350   -.577    .862   -.883
+  117SOL    HW2  351   -.465    .770   -.808
+  118SOL     OW  352   -.653    .525    .275
+  118SOL    HW1  353   -.640    .441    .329
+  118SOL    HW2  354   -.682    .599    .335
+  119SOL     OW  355    .307    .213   -.631
+  119SOL    HW1  356    .284    .250   -.541
+  119SOL    HW2  357    .277    .118   -.637
+  120SOL     OW  358    .037   -.552   -.580
+  120SOL    HW1  359    .090   -.601   -.512
+  120SOL    HW2  360    .059   -.454   -.575
+  121SOL     OW  361    .732    .634   -.798
+  121SOL    HW1  362    .791    .608   -.874
+  121SOL    HW2  363    .704    .730   -.809
+  122SOL     OW  364   -.134   -.927   -.008
+  122SOL    HW1  365   -.180   -.934   -.097
+  122SOL    HW2  366   -.196   -.883    .058
+  123SOL     OW  367    .307    .063    .618
+  123SOL    HW1  368    .296    .157    .651
+  123SOL    HW2  369    .302   -.000    .695
+  124SOL     OW  370   -.240    .367    .374
+  124SOL    HW1  371   -.238    .291    .438
+  124SOL    HW2  372   -.288    .444    .414
+  125SOL     OW  373   -.839    .766   -.896
+  125SOL    HW1  374   -.824    .787   -.800
+  125SOL    HW2  375   -.869    .671   -.905
+  126SOL     OW  376   -.882   -.289   -.162
+  126SOL    HW1  377   -.902   -.245   -.250
+  126SOL    HW2  378   -.843   -.380   -.178
+  127SOL     OW  379   -.003   -.344   -.257
+  127SOL    HW1  380    .011   -.317   -.352
+  127SOL    HW2  381    .080   -.322   -.204
+  128SOL     OW  382    .350    .898   -.058
+  128SOL    HW1  383    .426    .942   -.010
+  128SOL    HW2  384    .385    .851   -.140
+  129SOL     OW  385   -.322    .274    .125
+  129SOL    HW1  386   -.383    .199    .148
+  129SOL    HW2  387   -.300    .326    .208
+  130SOL     OW  388   -.559    .838    .042
+  130SOL    HW1  389   -.525    .745    .057
+  130SOL    HW2  390   -.541    .865   -.053
+  131SOL     OW  391   -.794   -.529    .849
+  131SOL    HW1  392   -.787   -.613    .794
+  131SOL    HW2  393   -.732   -.460    .813
+  132SOL     OW  394    .319    .810   -.913
+  132SOL    HW1  395    .412    .846   -.908
+  132SOL    HW2  396    .313    .725   -.861
+  133SOL     OW  397    .339    .509   -.856
+  133SOL    HW1  398    .287    .426   -.873
+  133SOL    HW2  399    .416    .514   -.920
+  134SOL     OW  400    .511    .415   -.054
+  134SOL    HW1  401    .493    .460    .034
+  134SOL    HW2  402    .553    .480   -.117
+  135SOL     OW  403   -.724    .380   -.184
+  135SOL    HW1  404   -.769    .443   -.120
+  135SOL    HW2  405   -.631    .411   -.201
+  136SOL     OW  406   -.702    .207   -.385
+  136SOL    HW1  407   -.702    .271   -.308
+  136SOL    HW2  408   -.674    .255   -.468
+  137SOL     OW  409    .008   -.536    .200
+  137SOL    HW1  410   -.085   -.515    .169
+  137SOL    HW2  411    .018   -.635    .213
+  138SOL     OW  412    .088   -.061    .927
+  138SOL    HW1  413    .046   -.147    .900
+  138SOL    HW2  414    .182   -.058    .893
+  139SOL     OW  415    .504   -.294    .910
+  139SOL    HW1  416    .570   -.220    .919
+  139SOL    HW2  417    .548   -.373    .868
+  140SOL     OW  418   -.860    .796   -.624
+  140SOL    HW1  419   -.819    .764   -.538
+  140SOL    HW2  420   -.956    .769   -.627
+  141SOL     OW  421    .040    .544   -.748
+  141SOL    HW1  422    .125    .511   -.789
+  141SOL    HW2  423    .053    .559   -.650
+  142SOL     OW  424    .189    .520   -.140
+  142SOL    HW1  425    .248    .480   -.210
+  142SOL    HW2  426    .131    .591   -.181
+  143SOL     OW  427   -.493   -.912   -.202
+  143SOL    HW1  428   -.454   -.823   -.182
+  143SOL    HW2  429   -.483   -.932   -.299
+  144SOL     OW  430    .815    .572    .325
+  144SOL    HW1  431    .822    .483    .279
+  144SOL    HW2  432    .721    .606    .317
+  145SOL     OW  433   -.205    .604   -.656
+  145SOL    HW1  434   -.243    .535   -.594
+  145SOL    HW2  435   -.123    .568   -.700
+  146SOL     OW  436    .252   -.298   -.118
+  146SOL    HW1  437    .222   -.241   -.042
+  146SOL    HW2  438    .245   -.395   -.092
+  147SOL     OW  439    .671    .464   -.593
+  147SOL    HW1  440    .637    .375   -.623
+  147SOL    HW2  441    .697    .518   -.673
+  148SOL     OW  442    .930   -.184   -.397
+  148SOL    HW1  443    .906   -.202   -.492
+  148SOL    HW2  444    .960   -.090   -.387
+  149SOL     OW  445    .473    .500    .191
+  149SOL    HW1  446    .534    .580    .195
+  149SOL    HW2  447    .378    .531    .198
+  150SOL     OW  448    .159   -.725   -.396
+  150SOL    HW1  449    .181   -.786   -.320
+  150SOL    HW2  450    .169   -.774   -.482
+  151SOL     OW  451   -.515   -.803   -.628
+  151SOL    HW1  452   -.491   -.866   -.702
+  151SOL    HW2  453   -.605   -.763   -.646
+  152SOL     OW  454   -.560    .855    .309
+  152SOL    HW1  455   -.646    .824    .351
+  152SOL    HW2  456   -.564    .841    .210
+  153SOL     OW  457   -.103   -.115   -.708
+  153SOL    HW1  458   -.042   -.085   -.781
+  153SOL    HW2  459   -.141   -.204   -.730
+  154SOL     OW  460   -.610   -.131   -.734
+  154SOL    HW1  461   -.526   -.126   -.788
+  154SOL    HW2  462   -.633   -.227   -.716
+  155SOL     OW  463    .083   -.604   -.840
+  155SOL    HW1  464    .078   -.605   -.740
+  155SOL    HW2  465    .000   -.645   -.878
+  156SOL     OW  466    .688   -.200   -.146
+  156SOL    HW1  467    .632   -.119   -.137
+  156SOL    HW2  468    .740   -.196   -.232
+  157SOL     OW  469    .903    .086    .133
+  157SOL    HW1  470    .954    .087    .047
+  157SOL    HW2  471    .959    .044    .204
+  158SOL     OW  472   -.136    .135    .523
+  158SOL    HW1  473   -.063    .118    .456
+  158SOL    HW2  474   -.167    .048    .561
+  159SOL     OW  475   -.474   -.289    .477
+  159SOL    HW1  476   -.407   -.277    .403
+  159SOL    HW2  477   -.514   -.200    .500
+  160SOL     OW  478    .130   -.068   -.011
+  160SOL    HW1  479    .089   -.142    .042
+  160SOL    HW2  480    .194   -.017    .047
+  161SOL     OW  481   -.582    .927    .672
+  161SOL    HW1  482   -.522    .846    .674
+  161SOL    HW2  483   -.542    .996    .612
+  162SOL     OW  484    .830   -.589   -.440
+  162SOL    HW1  485    .825   -.556   -.345
+  162SOL    HW2  486    .744   -.570   -.486
+  163SOL     OW  487    .672   -.246    .154
+  163SOL    HW1  488    .681   -.236    .055
+  163SOL    HW2  489    .632   -.335    .175
+  164SOL     OW  490   -.212   -.142   -.468
+  164SOL    HW1  491   -.159   -.132   -.552
+  164SOL    HW2  492   -.239   -.052   -.434
+  165SOL     OW  493   -.021    .175   -.899
+  165SOL    HW1  494    .018    .090   -.935
+  165SOL    HW2  495   -.119    .177   -.918
+  166SOL     OW  496    .263    .326    .720
+  166SOL    HW1  497    .184    .377    .686
+  166SOL    HW2  498    .254    .311    .818
+  167SOL     OW  499   -.668   -.250    .031
+  167SOL    HW1  500   -.662   -.343    .068
+  167SOL    HW2  501   -.727   -.250   -.049
+  168SOL     OW  502    .822   -.860   -.490
+  168SOL    HW1  503    .862   -.861   -.582
+  168SOL    HW2  504    .832   -.768   -.450
+  169SOL     OW  505    .916    .910    .291
+  169SOL    HW1  506    .979    .948    .223
+  169SOL    HW2  507    .956    .827    .330
+  170SOL     OW  508   -.358   -.255    .044
+  170SOL    HW1  509   -.450   -.218    .051
+  170SOL    HW2  510   -.320   -.235   -.046
+  171SOL     OW  511    .372   -.574   -.372
+  171SOL    HW1  512    .359   -.481   -.406
+  171SOL    HW2  513    .288   -.626   -.385
+  172SOL     OW  514   -.248   -.570   -.573
+  172SOL    HW1  515   -.188   -.567   -.493
+  172SOL    HW2  516   -.323   -.506   -.560
+  173SOL     OW  517   -.823   -.764    .696
+  173SOL    HW1  518   -.893   -.811    .750
+  173SOL    HW2  519   -.764   -.832    .653
+  174SOL     OW  520   -.848    .236   -.891
+  174SOL    HW1  521   -.856    .200   -.984
+  174SOL    HW2  522   -.850    .160   -.826
+  175SOL     OW  523    .590   -.375    .491
+  175SOL    HW1  524    .632   -.433    .421
+  175SOL    HW2  525    .546   -.296    .447
+  176SOL     OW  526   -.153    .385   -.481
+  176SOL    HW1  527   -.080    .454   -.477
+  176SOL    HW2  528   -.125    .310   -.540
+  177SOL     OW  529    .255   -.514    .290
+  177SOL    HW1  530    .159   -.513    .263
+  177SOL    HW2  531    .267   -.461    .374
+  178SOL     OW  532    .105   -.849   -.136
+  178SOL    HW1  533    .028   -.882   -.082
+  178SOL    HW2  534    .190   -.879   -.094
+  179SOL     OW  535    .672    .203   -.373
+  179SOL    HW1  536    .762    .187   -.413
+  179SOL    HW2  537    .680    .208   -.274
+  180SOL     OW  538    .075    .345    .033
+  180SOL    HW1  539   -.017    .317    .004
+  180SOL    HW2  540    .106    .422   -.023
+  181SOL     OW  541   -.422    .856   -.464
+  181SOL    HW1  542   -.479    .908   -.527
+  181SOL    HW2  543   -.326    .868   -.488
+  182SOL     OW  544    .072    .166    .318
+  182SOL    HW1  545    .055    .249    .264
+  182SOL    HW2  546    .162    .129    .296
+  183SOL     OW  547   -.679   -.527    .119
+  183SOL    HW1  548   -.778   -.538    .121
+  183SOL    HW2  549   -.645   -.512    .212
+  184SOL     OW  550    .613    .842   -.431
+  184SOL    HW1  551    .669    .923   -.448
+  184SOL    HW2  552    .672    .762   -.428
+  185SOL     OW  553   -.369   -.095   -.903
+  185SOL    HW1  554   -.336   -.031   -.972
+  185SOL    HW2  555   -.303   -.101   -.828
+  186SOL     OW  556    .716    .565   -.154
+  186SOL    HW1  557    .735    .630   -.080
+  186SOL    HW2  558    .776    .485   -.145
+  187SOL     OW  559   -.412   -.642   -.229
+  187SOL    HW1  560   -.421   -.652   -.130
+  187SOL    HW2  561   -.316   -.649   -.255
+  188SOL     OW  562    .390   -.121   -.302
+  188SOL    HW1  563    .299   -.080   -.304
+  188SOL    HW2  564    .383   -.215   -.270
+  189SOL     OW  565   -.188    .883   -.608
+  189SOL    HW1  566   -.215    .794   -.645
+  189SOL    HW2  567   -.187    .951   -.681
+  190SOL     OW  568   -.637    .325    .449
+  190SOL    HW1  569   -.572    .251    .438
+  190SOL    HW2  570   -.617    .375    .533
+  191SOL     OW  571    .594    .745    .652
+  191SOL    HW1  572    .644    .830    .633
+  191SOL    HW2  573    .506    .747    .604
+  192SOL     OW  574   -.085    .342   -.220
+  192SOL    HW1  575   -.102    .373   -.314
+  192SOL    HW2  576   -.169    .305   -.182
+  193SOL     OW  577   -.132   -.928   -.345
+  193SOL    HW1  578   -.094   -.837   -.330
+  193SOL    HW2  579   -.140   -.945   -.444
+  194SOL     OW  580    .859   -.488    .016
+  194SOL    HW1  581    .813   -.473    .104
+  194SOL    HW2  582    .903   -.403   -.014
+  195SOL     OW  583    .661   -.072   -.909
+  195SOL    HW1  584    .615    .016   -.922
+  195SOL    HW2  585    .760   -.060   -.916
+  196SOL     OW  586   -.454   -.011   -.142
+  196SOL    HW1  587   -.550   -.022   -.169
+  196SOL    HW2  588   -.398   -.078   -.190
+  197SOL     OW  589    .859   -.906    .861
+  197SOL    HW1  590    .913   -.975    .909
+  197SOL    HW2  591    .827   -.837    .927
+  198SOL     OW  592   -.779   -.878    .087
+  198SOL    HW1  593   -.802   -.825    .005
+  198SOL    HW2  594   -.698   -.934    .068
+  199SOL     OW  595   -.001   -.293    .851
+  199SOL    HW1  596   -.072   -.305    .781
+  199SOL    HW2  597    .000   -.372    .911
+  200SOL     OW  598    .221   -.548   -.018
+  200SOL    HW1  599    .156   -.621   -.039
+  200SOL    HW2  600    .225   -.534    .080
+  201SOL     OW  601    .079   -.622    .653
+  201SOL    HW1  602    .078   -.669    .741
+  201SOL    HW2  603    .161   -.650    .602
+  202SOL     OW  604    .672   -.471   -.238
+  202SOL    HW1  605    .594   -.521   -.200
+  202SOL    HW2  606    .669   -.376   -.207
+  203SOL     OW  607   -.038    .192   -.635
+  203SOL    HW1  608   -.042    .102   -.591
+  203SOL    HW2  609   -.035    .181   -.734
+  204SOL     OW  610    .428    .424    .520
+  204SOL    HW1  611    .458    .352    .458
+  204SOL    HW2  612    .389    .384    .603
+  205SOL     OW  613   -.157   -.375   -.758
+  205SOL    HW1  614   -.250   -.400   -.785
+  205SOL    HW2  615   -.131   -.425   -.676
+  206SOL     OW  616    .317    .547   -.582
+  206SOL    HW1  617    .355    .488   -.510
+  206SOL    HW2  618    .357    .521   -.670
+  207SOL     OW  619    .812   -.276    .687
+  207SOL    HW1  620    .844   -.266    .593
+  207SOL    HW2  621    .733   -.338    .689
+  208SOL     OW  622   -.438    .214   -.750
+  208SOL    HW1  623   -.386    .149   -.695
+  208SOL    HW2  624   -.487    .277   -.689
+  209SOL     OW  625   -.861    .034   -.708
+  209SOL    HW1  626   -.924   -.038   -.739
+  209SOL    HW2  627   -.768   -.002   -.708
+  210SOL     OW  628    .770   -.532    .301
+  210SOL    HW1  629    .724   -.619    .318
+  210SOL    HW2  630    .861   -.535    .342
+  211SOL     OW  631    .618   -.295   -.578
+  211SOL    HW1  632    .613   -.213   -.521
+  211SOL    HW2  633    .707   -.298   -.623
+  212SOL     OW  634   -.510    .052    .168
+  212SOL    HW1  635   -.475    .011    .084
+  212SOL    HW2  636   -.600    .014    .188
+  213SOL     OW  637   -.562    .453    .691
+  213SOL    HW1  638   -.621    .533    .695
+  213SOL    HW2  639   -.547    .418    .784
+  214SOL     OW  640   -.269    .221    .882
+  214SOL    HW1  641   -.353    .220    .936
+  214SOL    HW2  642   -.267    .304    .826
+  215SOL     OW  643    .039   -.785    .300
+  215SOL    HW1  644    .138   -.796    .291
+  215SOL    HW2  645   -.001   -.871    .332
+  216SOL     OW  646    .875   -.216    .337
+  216SOL    HW1  647    .798   -.251    .283
+  216SOL    HW2  648    .843   -.145    .399
+   1.86206   1.86206   1.86206
index 04a490543825cb794406527dfba300134377ae80..21579beb2aa3dfbba310891368541270799706b8 100644 (file)
@@ -105,7 +105,6 @@ set(SYMLINK_NAMES
     g_potential
     g_principal
     g_rama
-    g_rdf
     g_rms
     g_rmsdist
     g_rmsf
index 720122cba70ab45d8a5adc810d0a9b18a2d087df..71a2394de3c8a16bc39695005cb0e47acc3a6e17 100644 (file)
@@ -363,8 +363,6 @@ void registerLegacyModules(gmx::CommandLineModuleManager *manager)
                    "Calculate principal axes of inertia for a group of atoms");
     registerModule(manager, &gmx_rama, "rama",
                    "Compute Ramachandran plots");
-    registerModule(manager, &gmx_rdf, "rdf",
-                   "Calculate radial distribution functions");
     registerModule(manager, &gmx_rms, "rms",
                    "Calculate RMSDs with a reference structure and RMSD matrices");
     registerModule(manager, &gmx_rmsdist, "rmsdist",