Add occupied fraction and PDB output to g_ana select.

- Add output option -of to write out the occupied fraction.
- Add output option -ofpdb to write the same information into a PDB
  file. Option -pdbatoms can be used to control which atoms to write
  out.
- Reorganize error checking for options that only work with a single
  selection or only make sense with dynamic selections.

Supporting changes:
- Add PDB file type to FileNameOption.
- Make write_trxframe_indexed() take the index as a const array.
- Add ConstArrayRef::data() following std::vector::data() from C++11.

The PDB output can also be used to more easily check what is selected by
a selection.

Change-Id: Icabab4b2c88cc21e4d5ebb672a1471c254acb05e

diff --git a/src/gromacs/gmxlib/confio.c b/src/gromacs/gmxlib/confio.c
index 886983a6d513dda269c8ad5b106ebaca6fd74ac5..0af17e2bb4d6919f3b5ef0d9e7e9355372f83b46 100644 (file)
--- a/src/gromacs/gmxlib/confio.c
+++ b/src/gromacs/gmxlib/confio.c
@@ -288,7 +288,7 @@ int read_g96_conf(FILE *fp,const char *infile,t_trxframe *fr, char *line)
 }
 
 void write_g96_conf(FILE *out,t_trxframe *fr,
-                   int nindex,atom_id *index)
+                    int nindex,const atom_id *index)
 {
   t_atoms *atoms;
   int nout,i,a;
@@ -993,7 +993,7 @@ static void write_hconf_box(FILE *out,int pr,matrix box)
 }
 
 void write_hconf_indexed_p(FILE *out,const char *title,t_atoms *atoms,
-                          int nx,atom_id index[], int pr,
+                          int nx,const atom_id index[], int pr,
                           rvec *x,rvec *v,matrix box)
 {
   char resnm[6],nm[6],format[100];

diff --git a/src/gromacs/gmxlib/pdbio.c b/src/gromacs/gmxlib/pdbio.c
index 4c7ea5b92766ee3819cd2230d1161c4fafdc7b3b..a7ead4ac993dd8d90f503ff078db2fe293437643 100644 (file)
--- a/src/gromacs/gmxlib/pdbio.c
+++ b/src/gromacs/gmxlib/pdbio.c
@@ -219,7 +219,7 @@ static void read_cryst1(char *line,int *ePBC,matrix box)
 void write_pdbfile_indexed(FILE *out,const char *title,
                           t_atoms *atoms,rvec x[],
                           int ePBC,matrix box,char chainid,
-                          int model_nr, atom_id nindex, atom_id index[],
+                          int model_nr, atom_id nindex, const atom_id index[],
                           gmx_conect conect, gmx_bool bTerSepChains)
 {
   gmx_conect_t *gc = (gmx_conect_t *)conect;

diff --git a/src/gromacs/gmxlib/trxio.c b/src/gromacs/gmxlib/trxio.c
index 14a1ff592d0c28bbfe8ebb935efdc8ea13dcc267..50d1ed01171df9a982b5962844c959abe34fe3fc 100644 (file)
--- a/src/gromacs/gmxlib/trxio.c
+++ b/src/gromacs/gmxlib/trxio.c
@@ -189,7 +189,7 @@ void set_trxframe_ePBC(t_trxframe *fr,int ePBC)
 }
 
 int write_trxframe_indexed(t_trxstatus *status,t_trxframe *fr,int nind,
-                           atom_id *ind, gmx_conect gc)
+                           const atom_id *ind, gmx_conect gc)
 {
   char title[STRLEN];
   rvec *xout=NULL,*vout=NULL,*fout=NULL;
@@ -354,7 +354,7 @@ int write_trxframe(t_trxstatus *status,t_trxframe *fr,gmx_conect gc)
   return 0;
 }
 
-int write_trx(t_trxstatus *status,int nind,atom_id *ind,t_atoms *atoms,
+int write_trx(t_trxstatus *status,int nind,const atom_id *ind,t_atoms *atoms,
              int step,real time,matrix box,rvec x[],rvec *v,
              gmx_conect gc)
 {

diff --git a/src/gromacs/legacyheaders/confio.h b/src/gromacs/legacyheaders/confio.h
index a8238babc66c603e79b27bec7659759fed92115c..2b2372545fada9d30c36ccc98f3cf1208dd477ce 100644 (file)
--- a/src/gromacs/legacyheaders/confio.h
+++ b/src/gromacs/legacyheaders/confio.h
@@ -60,7 +60,7 @@ int read_g96_conf(FILE *fp,const char *infile,t_trxframe *fr, char *line);
  * title, atoms, x, v can all be NULL, in which case they won't be read *
  * line holds the previous line for trajectory reading                  */
 
-void write_g96_conf(FILE *out,t_trxframe *fr,int nindex,atom_id *index);
+void write_g96_conf(FILE *out,t_trxframe *fr,int nindex,const atom_id *index);
 /* write a Gromos96 coordinate file or trajectory frame *
  * index can be NULL                                    */
 
@@ -69,9 +69,9 @@ int gro_first_x_or_v(FILE *status,t_trxframe *fr);
 /* read first/next x and/or v frame from gro file */
 
 void write_hconf_indexed_p(FILE *out,const char *title,t_atoms *atoms,
-                                 int nx,atom_id index[],int ndec,
-                                 rvec *x,rvec *v,matrix box);
-               
+                           int nx,const atom_id index[],int ndec,
+                           rvec *x,rvec *v,matrix box);
+
 void write_hconf_p(FILE *out,const char *title,t_atoms *atoms, int ndec,
                          rvec *x,rvec *v,matrix box); 
 /* Write a Gromos file with precision ndec: number of decimal places in x,

diff --git a/src/gromacs/legacyheaders/pdbio.h b/src/gromacs/legacyheaders/pdbio.h
index f1dfd8fb9093aa341e241b86c9f70814a534ee7a..5230688f1e1a621dca5d1f9d577a6e264fdef173 100644 (file)
--- a/src/gromacs/legacyheaders/pdbio.h
+++ b/src/gromacs/legacyheaders/pdbio.h
@@ -75,9 +75,9 @@ void gmx_write_pdb_box(FILE *out,int ePBC,matrix box);
  */
 
 void write_pdbfile_indexed(FILE *out,const char *title,t_atoms *atoms,
-                                 rvec x[],int ePBC,matrix box,char chain,
-                                 int model_nr,atom_id nindex,atom_id index[],
-                                 gmx_conect conect,gmx_bool bTerSepChains);
+                           rvec x[],int ePBC,matrix box,char chain,
+                           int model_nr,atom_id nindex,const atom_id index[],
+                           gmx_conect conect,gmx_bool bTerSepChains);
 /* REALLY low level */
 
 void write_pdbfile(FILE *out,const char *title,t_atoms *atoms,

diff --git a/src/gromacs/legacyheaders/statutil.h b/src/gromacs/legacyheaders/statutil.h
index 9620f2c28733b1509d43726edc16446bbff771ea..24701c0c19158d907220aa2e2b395bad3ff175c6 100644 (file)
--- a/src/gromacs/legacyheaders/statutil.h
+++ b/src/gromacs/legacyheaders/statutil.h
@@ -111,7 +111,7 @@ int nframes_read(t_trxstatus *status);
 /* Returns the number of frames read from the trajectory */
 
 int write_trxframe_indexed(t_trxstatus *status,t_trxframe *fr,int nind,
-                           atom_id *ind, gmx_conect gc);
+                           const atom_id *ind, gmx_conect gc);
 /* Write an indexed frame to a TRX file, see write_trxframe. gc may be NULL */
 
 int write_trxframe(t_trxstatus *status,t_trxframe *fr,gmx_conect gc);
@@ -124,7 +124,7 @@ int write_trxframe(t_trxstatus *status,t_trxframe *fr,gmx_conect gc);
  * gc is important for pdb file writing only and may be NULL.
  */
 
-int write_trx(t_trxstatus *status,int nind,atom_id *ind,t_atoms *atoms,
+int write_trx(t_trxstatus *status,int nind,const atom_id *ind,t_atoms *atoms,
               int step,real time,matrix box,rvec x[],rvec *v,
               gmx_conect gc);
 /* Write an indexed frame to a TRX file.

diff --git a/src/gromacs/options/filenameoption.cpp b/src/gromacs/options/filenameoption.cpp
index cee1e04b6344e19401a0287ea8f69d18036f0747..798c1887d9e302f36d257022e4d71f3f196aade1 100644 (file)
--- a/src/gromacs/options/filenameoption.cpp
+++ b/src/gromacs/options/filenameoption.cpp
@@ -184,6 +184,7 @@ FileTypeRegistry::FileTypeRegistry()
     };
     registerType(eftTopology,    topExtensions);
     registerType(eftTrajectory,  trajExtensions);
+    registerType(eftPDB,         ".pdb");
     registerType(eftIndex,       ".ndx");
     registerType(eftPlot,        ".xvg");
     registerType(eftGenericData, ".dat");

diff --git a/src/gromacs/options/optionfiletype.h b/src/gromacs/options/optionfiletype.h
index 675cf57fef67de258846308ead8ba75b444dc1b6..f263eb548331d1de33edd1a0992c8a3e4305390a 100644 (file)
--- a/src/gromacs/options/optionfiletype.h
+++ b/src/gromacs/options/optionfiletype.h
@@ -51,6 +51,7 @@ enum OptionFileType {
     eftUnknown,
     eftTopology,
     eftTrajectory,
+    eftPDB,
     eftIndex,
     eftPlot,
     eftGenericData,

diff --git a/src/gromacs/trajectoryanalysis/modules/select.cpp b/src/gromacs/trajectoryanalysis/modules/select.cpp
index 76c8d300803d034f6ac402591f975c1aac9b7a36..93a686524f7df093ee20e2d69f4442bfb28e5d50 100644 (file)
--- a/src/gromacs/trajectoryanalysis/modules/select.cpp
+++ b/src/gromacs/trajectoryanalysis/modules/select.cpp
@@ -37,20 +37,21 @@
  */
 #include "select.h"
 
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include <cstdio>
+#include <cstring>
 
 #include <algorithm>
-#include <cstdio>
 #include <string>
 #include <vector>
 
-#include "gmxfio.h"
+#include "gromacs/legacyheaders/gmxfio.h"
+#include "gromacs/legacyheaders/smalloc.h"
+#include "gromacs/legacyheaders/statutil.h"
 
 #include "gromacs/analysisdata/analysisdata.h"
 #include "gromacs/analysisdata/dataframe.h"
 #include "gromacs/analysisdata/datamodule.h"
+#include "gromacs/analysisdata/modules/average.h"
 #include "gromacs/analysisdata/modules/plot.h"
 #include "gromacs/options/basicoptions.h"
 #include "gromacs/options/filenameoption.h"
@@ -59,6 +60,7 @@
 #include "gromacs/selection/selectionoption.h"
 #include "gromacs/trajectoryanalysis/analysissettings.h"
 #include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/stringutil.h"
 
 namespace gmx
@@ -245,8 +247,9 @@ const char Select::shortDescription[] =
 
 Select::Select()
     : TrajectoryAnalysisModule(name, shortDescription),
+      selOpt_(NULL),
       bDump_(false), bTotNorm_(false), bFracNorm_(false), bResInd_(false),
-      top_(NULL)
+      top_(NULL), occupancyModule_(new AnalysisDataAverageModule())
 {
     registerAnalysisDataset(&sdata_, "size");
     registerAnalysisDataset(&cdata_, "cfrac");
@@ -254,6 +257,8 @@ Select::Select()
     idata_.setMultipoint(true);
     registerAnalysisDataset(&idata_, "index");
     registerAnalysisDataset(&mdata_, "mask");
+    occupancyModule_->setXAxis(1.0, 1.0);
+    registerBasicDataset(occupancyModule_.get(), "occupancy");
 }
 
 
@@ -310,7 +315,21 @@ Select::initOptions(Options *options, TrajectoryAnalysisSettings * /*settings*/)
         "one frame, and contains either 0/1 for each atom/residue/molecule",
         "possibly selected. 1 stands for the atom/residue/molecule being",
         "selected for the current frame, 0 for not selected.",
-        "With [TT]-dump[tt], the frame time is omitted from the output."
+        "With [TT]-dump[tt], the frame time is omitted from the output.[PAR]",
+        "With [TT]-of[tt], the occupancy fraction of each position (i.e.,",
+        "the fraction of frames where the position is selected) is",
+        "printed.[PAR]",
+        "With [TT]-ofpdb[tt], a PDB file is written out where the occupancy",
+        "column is filled with the occupancy fraction of each atom in the",
+        "selection. The coordinates in the PDB file will be those from the",
+        "input topology. [TT]-pdbatoms[tt] can be used to control which atoms",
+        "appear in the output PDB file: with [TT]all[tt] all atoms are",
+        "present, with [TT]maxsel[tt] all atoms possibly selected by the",
+        "selection are present, and with [TT]selected[tt] only atoms that are",
+        "selected at least in one frame are present.[PAR]",
+        "With [TT]-om[tt], [TT]-of[tt] and [TT]-ofpdb[tt], only one selection",
+        "can be provided. [TT]-om[tt] and [TT]-of[tt] only accept dynamic",
+        "selections."
     };
 
     options->setDescription(concatenateStrings(desc));
@@ -319,7 +338,7 @@ Select::initOptions(Options *options, TrajectoryAnalysisSettings * /*settings*/)
                            .store(&fnSize_).defaultBasename("size")
                            .description("Number of positions in each selection"));
     options->addOption(FileNameOption("oc").filetype(eftPlot).outputFile()
-                           .store(&fnFrac_).defaultBasename("frac")
+                           .store(&fnFrac_).defaultBasename("cfrac")
                            .description("Covered fraction for each selection"));
     options->addOption(FileNameOption("oi").filetype(eftGenericData).outputFile()
                            .store(&fnIndex_).defaultBasename("index")
@@ -330,8 +349,14 @@ Select::initOptions(Options *options, TrajectoryAnalysisSettings * /*settings*/)
     options->addOption(FileNameOption("om").filetype(eftPlot).outputFile()
                            .store(&fnMask_).defaultBasename("mask")
                            .description("Mask for selected positions"));
-
-    options->addOption(SelectionOption("select").storeVector(&sel_)
+    options->addOption(FileNameOption("of").filetype(eftPlot).outputFile()
+                           .store(&fnOccupancy_).defaultBasename("occupancy")
+                           .description("Occupied fraction for selected positions"));
+    options->addOption(FileNameOption("ofpdb").filetype(eftPDB).outputFile()
+                           .store(&fnPDB_).defaultBasename("occupancy")
+                           .description("PDB file with occupied fraction for selected positions"));
+
+    selOpt_ = options->addOption(SelectionOption("select").storeVector(&sel_)
         .required().multiValue()
         .description("Selections to analyze"));
 
@@ -344,17 +369,37 @@ Select::initOptions(Options *options, TrajectoryAnalysisSettings * /*settings*/)
     const char *const cResNumberEnum[] = { "number", "index", NULL };
     options->addOption(StringOption("resnr").store(&resNumberType_)
         .enumValue(cResNumberEnum).defaultEnumIndex(0)
-        .description("Residue number output type"));
+        .description("Residue number output type with -oi and -on"));
+    const char *const cPDBAtomsEnum[] = { "all", "maxsel", "selected", NULL };
+    options->addOption(StringOption("pdbatoms").store(&pdbAtoms_)
+        .enumValue(cPDBAtomsEnum).defaultEnumIndex(0)
+        .description("Atoms to write with -ofpdb"));
 }
 
+void
+Select::optionsFinished(Options * /*options*/,
+                        TrajectoryAnalysisSettings *settings)
+{
+    if (!fnPDB_.empty())
+    {
+        settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
+        settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
+    }
+    if ((!fnIndex_.empty() && bDump_)
+        || !fnMask_.empty() || !fnOccupancy_.empty() || !fnPDB_.empty())
+    {
+        selOpt_->setValueCount(1);
+    }
+}
 
 void
 Select::initAnalysis(const TrajectoryAnalysisSettings &settings,
                      const TopologyInformation &top)
 {
-    if (!fnIndex_.empty() && bDump_ && sel_.size() > 1U)
+    if (!sel_[0].isDynamic() && (!fnMask_.empty() || !fnOccupancy_.empty()))
     {
-        GMX_THROW(InconsistentInputError("With -oi and -dump, there can be only one selection"));
+        GMX_THROW(InconsistentInputError(
+                    "-om or -of are not meaningful with a static selection"));
     }
     bResInd_ = (resNumberType_ == "index");
 
@@ -424,32 +469,32 @@ Select::initAnalysis(const TrajectoryAnalysisSettings &settings,
     }
 
     mdata_.setColumnCount(sel_[0].posCount());
+    mdata_.addModule(occupancyModule_);
     if (!fnMask_.empty())
     {
-        if (sel_.size() > 1U)
-        {
-            fprintf(stderr, "WARNING: the mask (-om) will only be written for the first group\n");
-        }
-        if (!sel_[0].isDynamic())
-        {
-            fprintf(stderr, "WARNING: will not write the mask (-om) for a static selection\n");
-        }
-        else
-        {
-            AnalysisDataPlotModulePointer plot(
-                new AnalysisDataPlotModule(settings.plotSettings()));
-            plot->setFileName(fnMask_);
-            plot->setPlainOutput(bDump_);
-            plot->setOmitX(bDump_);
-            plot->setTitle("Selection mask");
-            plot->setXAxisIsTime();
-            plot->setYLabel("Occupancy");
-            plot->setYFormat(1, 0);
-            mdata_.addModule(plot);
-        }
+        AnalysisDataPlotModulePointer plot(
+            new AnalysisDataPlotModule(settings.plotSettings()));
+        plot->setFileName(fnMask_);
+        plot->setPlainOutput(bDump_);
+        plot->setOmitX(bDump_);
+        plot->setTitle("Selection mask");
+        plot->setXAxisIsTime();
+        plot->setYLabel("Occupancy");
+        plot->setYFormat(1, 0);
+        mdata_.addModule(plot);
+    }
+    if (!fnOccupancy_.empty())
+    {
+        AnalysisDataPlotModulePointer plot(
+            new AnalysisDataPlotModule(settings.plotSettings()));
+        plot->setFileName(fnOccupancy_);
+        plot->setTitle("Fraction of time selection matches");
+        plot->setXLabel("Selected position");
+        plot->setYLabel("Occupied fraction");
+        occupancyModule_->addColumnModule(0, 1, plot);
     }
 
-    top_ = top.topology();
+    top_ = &top;
 }
 
 
@@ -462,6 +507,7 @@ Select::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
     AnalysisDataHandle idh = pdata->dataHandle(idata_);
     AnalysisDataHandle mdh = pdata->dataHandle(mdata_);
     const SelectionList &sel = pdata->parallelSelections(sel_);
+    t_topology *top = top_->topology();
 
     sdh.startFrame(frnr, fr.time);
     for (size_t g = 0; g < sel.size(); ++g)
@@ -492,7 +538,7 @@ Select::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
             const SelectionPosition &p = sel[g].position(i);
             if (sel[g].type() == INDEX_RES && !bResInd_)
             {
-                idh.setPoint(1, top_->atoms.resinfo[p.mappedId()].nr);
+                idh.setPoint(1, top->atoms.resinfo[p.mappedId()].nr);
             }
             else
             {
@@ -525,6 +571,81 @@ Select::finishAnalysis(int /*nframes*/)
 void
 Select::writeOutput()
 {
+    if (!fnPDB_.empty())
+    {
+        GMX_RELEASE_ASSERT(top_->hasTopology(),
+                "Topology should have been loaded or an error given earlier");
+        t_atoms atoms;
+        atoms = top_->topology()->atoms;
+        t_pdbinfo *pdbinfo;
+        snew(pdbinfo, atoms.nr);
+        scoped_ptr_sfree pdbinfoGuard(pdbinfo);
+        if (atoms.pdbinfo != NULL)
+        {
+            std::memcpy(pdbinfo, atoms.pdbinfo, atoms.nr*sizeof(*pdbinfo));
+        }
+        atoms.pdbinfo = pdbinfo;
+        for (int i = 0; i < atoms.nr; ++i)
+        {
+            pdbinfo[i].occup = 0.0;
+        }
+        for (int i = 0; i < sel_[0].posCount(); ++i)
+        {
+            ConstArrayRef<int> atomIndices = sel_[0].position(i).atomIndices();
+            ConstArrayRef<int>::const_iterator ai;
+            for (ai = atomIndices.begin(); ai != atomIndices.end(); ++ai)
+            {
+                pdbinfo[*ai].occup = occupancyModule_->average(i);
+            }
+        }
+
+        t_trxframe fr;
+        clear_trxframe(&fr, TRUE);
+        fr.bAtoms = TRUE;
+        fr.atoms  = &atoms;
+        fr.bX     = TRUE;
+        fr.bBox   = TRUE;
+        top_->getTopologyConf(&fr.x, fr.box);
+
+        if (pdbAtoms_ == "all")
+        {
+            t_trxstatus *status = open_trx(fnPDB_.c_str(), "w");
+            write_trxframe(status, &fr, NULL);
+            close_trx(status);
+        }
+        else if (pdbAtoms_ == "maxsel")
+        {
+            ConstArrayRef<int> atomIndices = sel_[0].atomIndices();
+            t_trxstatus *status = open_trx(fnPDB_.c_str(), "w");
+            write_trxframe_indexed(status, &fr, atomIndices.size(),
+                                   atomIndices.data(), NULL);
+            close_trx(status);
+        }
+        else if (pdbAtoms_ == "selected")
+        {
+            std::vector<int> indices;
+            for (int i = 0; i < sel_[0].posCount(); ++i)
+            {
+                if (occupancyModule_->average(i) > 0)
+                {
+                    ConstArrayRef<int> atomIndices = sel_[0].position(i).atomIndices();
+                    ConstArrayRef<int>::const_iterator ai;
+                    for (ai = atomIndices.begin(); ai != atomIndices.end(); ++ai)
+                    {
+                        indices.push_back(*ai);
+                    }
+                }
+            }
+            t_trxstatus *status = open_trx(fnPDB_.c_str(), "w");
+            write_trxframe_indexed(status, &fr, indices.size(), &indices[0], NULL);
+            close_trx(status);
+        }
+        else
+        {
+            GMX_RELEASE_ASSERT(false,
+                    "Mismatch between -pdbatoms enum values and implementation");
+        }
+    }
 }
 
 } // namespace analysismodules

diff --git a/src/gromacs/trajectoryanalysis/modules/select.h b/src/gromacs/trajectoryanalysis/modules/select.h
index d4ca4e2f9470580c87abbb49056be3f1a22d7697..b0434017197aa63000e85fc6b23033f63056a4cd 100644 (file)
--- a/src/gromacs/trajectoryanalysis/modules/select.h
+++ b/src/gromacs/trajectoryanalysis/modules/select.h
@@ -43,11 +43,14 @@
 
 #include "../analysismodule.h"
 #include "gromacs/analysisdata/analysisdata.h"
+#include "gromacs/analysisdata/modules/average.h"
 #include "gromacs/selection/selection.h"
 
 namespace gmx
 {
 
+class SelectionOptionInfo;
+
 namespace analysismodules
 {
 
@@ -62,6 +65,8 @@ class Select : public TrajectoryAnalysisModule
 
         virtual void initOptions(Options *options,
                                  TrajectoryAnalysisSettings *settings);
+        virtual void optionsFinished(Options *options,
+                                     TrajectoryAnalysisSettings *settings);
         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
                                   const TopologyInformation &top);
 
@@ -73,24 +78,29 @@ class Select : public TrajectoryAnalysisModule
 
     private:
         SelectionList            sel_;
+        SelectionOptionInfo     *selOpt_;
 
         std::string              fnSize_;
         std::string              fnFrac_;
         std::string              fnIndex_;
         std::string              fnNdx_;
         std::string              fnMask_;
+        std::string              fnOccupancy_;
+        std::string              fnPDB_;
         bool                     bDump_;
         bool                     bTotNorm_;
         bool                     bFracNorm_;
         bool                     bResInd_;
         std::string              resNumberType_;
+        std::string              pdbAtoms_;
 
-        t_topology              *top_;
+        const TopologyInformation *top_;
         std::vector<int>         totsize_;
         AnalysisData             sdata_;
         AnalysisData             cdata_;
         AnalysisData             idata_;
         AnalysisData             mdata_;
+        AnalysisDataAverageModulePointer occupancyModule_;
 };
 
 } // namespace analysismodules

diff --git a/src/gromacs/trajectoryanalysis/tests/refdata/SelectModuleTest_BasicTest.xml b/src/gromacs/trajectoryanalysis/tests/refdata/SelectModuleTest_BasicTest.xml
index 8349ca5e3f6a57c1036b399610aabc9f4f968587..eef1532f568da7549f2c4e35d787ead220bcd215 100644 (file)
--- a/src/gromacs/trajectoryanalysis/tests/refdata/SelectModuleTest_BasicTest.xml
+++ b/src/gromacs/trajectoryanalysis/tests/refdata/SelectModuleTest_BasicTest.xml
@@ -135,6 +135,137 @@
           </DataValue>
         </Sequence>
       </DataFrame>
+      <DataFrame Name="Frame1">
+        <Real Name="X">0.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">1</Int>
+          <Int Name="FirstColumn">0</Int>
+          <Int Name="LastColumn">0</Int>
+          <DataValue>
+            <Real Name="Value">8.000000</Real>
+          </DataValue>
+        </Sequence>
+        <Sequence Name="Y">
+          <Int Name="Length">1</Int>
+          <Int Name="FirstColumn">1</Int>
+          <Int Name="LastColumn">1</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+        </Sequence>
+        <Sequence Name="Y">
+          <Int Name="Length">1</Int>
+          <Int Name="FirstColumn">1</Int>
+          <Int Name="LastColumn">1</Int>
+          <DataValue>
+            <Real Name="Value">3.000000</Real>
+          </DataValue>
+        </Sequence>
+        <Sequence Name="Y">
+          <Int Name="Length">1</Int>
+          <Int Name="FirstColumn">1</Int>
+          <Int Name="LastColumn">1</Int>
+          <DataValue>
+            <Real Name="Value">4.000000</Real>
+          </DataValue>
+        </Sequence>
+        <Sequence Name="Y">
+          <Int Name="Length">1</Int>
+          <Int Name="FirstColumn">1</Int>
+          <Int Name="LastColumn">1</Int>
+          <DataValue>
+            <Real Name="Value">7.000000</Real>
+          </DataValue>
+        </Sequence>
+        <Sequence Name="Y">
+          <Int Name="Length">1</Int>
+          <Int Name="FirstColumn">1</Int>
+          <Int Name="LastColumn">1</Int>
+          <DataValue>
+            <Real Name="Value">8.000000</Real>
+          </DataValue>
+        </Sequence>
+        <Sequence Name="Y">
+          <Int Name="Length">1</Int>
+          <Int Name="FirstColumn">1</Int>
+          <Int Name="LastColumn">1</Int>
+          <DataValue>
+            <Real Name="Value">9.000000</Real>
+          </DataValue>
+        </Sequence>
+        <Sequence Name="Y">
+          <Int Name="Length">1</Int>
+          <Int Name="FirstColumn">1</Int>
+          <Int Name="LastColumn">1</Int>
+          <DataValue>
+            <Real Name="Value">12.000000</Real>
+          </DataValue>
+        </Sequence>
+        <Sequence Name="Y">
+          <Int Name="Length">1</Int>
+          <Int Name="FirstColumn">1</Int>
+          <Int Name="LastColumn">1</Int>
+          <DataValue>
+            <Real Name="Value">13.000000</Real>
+          </DataValue>
+        </Sequence>
+        <Sequence Name="Y">
+          <Int Name="Length">1</Int>
+          <Int Name="FirstColumn">0</Int>
+          <Int Name="LastColumn">0</Int>
+          <DataValue>
+            <Real Name="Value">6.000000</Real>
+          </DataValue>
+        </Sequence>
+        <Sequence Name="Y">
+          <Int Name="Length">1</Int>
+          <Int Name="FirstColumn">1</Int>
+          <Int Name="LastColumn">1</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+        </Sequence>
+        <Sequence Name="Y">
+          <Int Name="Length">1</Int>
+          <Int Name="FirstColumn">1</Int>
+          <Int Name="LastColumn">1</Int>
+          <DataValue>
+            <Real Name="Value">2.000000</Real>
+          </DataValue>
+        </Sequence>
+        <Sequence Name="Y">
+          <Int Name="Length">1</Int>
+          <Int Name="FirstColumn">1</Int>
+          <Int Name="LastColumn">1</Int>
+          <DataValue>
+            <Real Name="Value">3.000000</Real>
+          </DataValue>
+        </Sequence>
+        <Sequence Name="Y">
+          <Int Name="Length">1</Int>
+          <Int Name="FirstColumn">1</Int>
+          <Int Name="LastColumn">1</Int>
+          <DataValue>
+            <Real Name="Value">7.000000</Real>
+          </DataValue>
+        </Sequence>
+        <Sequence Name="Y">
+          <Int Name="Length">1</Int>
+          <Int Name="FirstColumn">1</Int>
+          <Int Name="LastColumn">1</Int>
+          <DataValue>
+            <Real Name="Value">8.000000</Real>
+          </DataValue>
+        </Sequence>
+        <Sequence Name="Y">
+          <Int Name="Length">1</Int>
+          <Int Name="FirstColumn">1</Int>
+          <Int Name="LastColumn">1</Int>
+          <DataValue>
+            <Real Name="Value">9.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
     </AnalysisData>
     <AnalysisData Name="mask">
       <DataFrame Name="Frame0">
@@ -188,6 +319,239 @@
           </DataValue>
         </Sequence>
       </DataFrame>
+      <DataFrame Name="Frame1">
+        <Real Name="X">0.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">15</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+    </AnalysisData>
+    <AnalysisData Name="occupancy">
+      <DataFrame Name="Frame0">
+        <Real Name="X">1.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame1">
+        <Real Name="X">2.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame2">
+        <Real Name="X">3.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame3">
+        <Real Name="X">4.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame4">
+        <Real Name="X">5.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame5">
+        <Real Name="X">6.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame6">
+        <Real Name="X">7.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame7">
+        <Real Name="X">8.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame8">
+        <Real Name="X">9.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame9">
+        <Real Name="X">10.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame10">
+        <Real Name="X">11.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame11">
+        <Real Name="X">12.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame12">
+        <Real Name="X">13.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame13">
+        <Real Name="X">14.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame14">
+        <Real Name="X">15.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
     </AnalysisData>
     <AnalysisData Name="size">
       <DataFrame Name="Frame0">
@@ -202,11 +566,24 @@
           </DataValue>
         </Sequence>
       </DataFrame>
+      <DataFrame Name="Frame1">
+        <Real Name="X">0.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">8.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">6.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
     </AnalysisData>
   </OutputData>
   <OutputFiles Name="Files">
     <String Name="-oi"><![CDATA[
       0.000    8    1    2    5    6    9   10   13   14    6    1    2    3    7    8    9
+      0.000    8    1    3    4    7    8    9   12   13    6    1    2    3    7    8    9
 ]]></String>
     <String Name="-on"><![CDATA[
 [ y_<_2.5_f0_t0.000 ]
@@ -214,6 +591,9 @@
 
 [ resname_RA ]
    1    2    3    7    8    9 
+
+[ y_<_2.5_f1_t0.000 ]
+   1    3    4    7    8    9   12   13 
 ]]></String>
   </OutputFiles>
 </ReferenceData>

diff --git a/src/gromacs/trajectoryanalysis/tests/refdata/SelectModuleTest_HandlesMaxPDBOutput.xml b/src/gromacs/trajectoryanalysis/tests/refdata/SelectModuleTest_HandlesMaxPDBOutput.xml
new file mode 100644 (file)
index 0000000..97ff96d
--- /dev/null
+++ b/src/gromacs/trajectoryanalysis/tests/refdata/SelectModuleTest_HandlesMaxPDBOutput.xml
@@ -0,0 +1,136 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <String Name="CommandLine">select -select 'resname RA RD and y &lt; 2.5' -pdbatoms maxsel</String>
+  <OutputData Name="Data">
+    <AnalysisData Name="occupancy">
+      <DataFrame Name="Frame0">
+        <Real Name="X">1.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame1">
+        <Real Name="X">2.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame2">
+        <Real Name="X">3.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame3">
+        <Real Name="X">4.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame4">
+        <Real Name="X">5.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame5">
+        <Real Name="X">6.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame6">
+        <Real Name="X">7.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame7">
+        <Real Name="X">8.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame8">
+        <Real Name="X">9.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+    </AnalysisData>
+  </OutputData>
+  <OutputFiles Name="Files">
+    <String Name="-ofpdb"><![CDATA[
+TITLE     frame t= 0.000
+REMARK    THIS IS A SIMULATION BOX
+CRYST1  100.000  100.000  100.000  90.00  90.00  90.00 P 1           1
+MODEL        1
+ATOM      1  CB   RA A   1      10.000  10.000   0.000  1.00  0.00            
+ATOM      2  S1   RA A   1      10.000  20.000   0.000  0.50  0.00            
+ATOM      3  S2   RA A   1      10.000  30.000   0.000  0.50  0.00            
+ATOM      7  CB   RA B   1      20.000  30.000   0.000  0.50  0.00            
+ATOM      8  S1   RA B   1      20.000  40.000   0.000  0.50  0.30            
+ATOM      9  S2   RA B   1      30.000  10.000   0.000  1.00  0.30            
+ATOM     13  CB   RD B   2      40.000  10.000   0.000  1.00  0.00            
+ATOM     14  S1   RD B   2      40.000  20.000   0.000  0.50  0.00            
+ATOM     15  S2   RD B   2      40.000  30.000   0.000  0.00  0.00            
+TER
+ENDMDL
+]]></String>
+  </OutputFiles>
+</ReferenceData>

diff --git a/src/gromacs/trajectoryanalysis/tests/refdata/SelectModuleTest_HandlesPDBOutputWithNonPDBInput.xml b/src/gromacs/trajectoryanalysis/tests/refdata/SelectModuleTest_HandlesPDBOutputWithNonPDBInput.xml
new file mode 100644 (file)
index 0000000..cb042f3
--- /dev/null
+++ b/src/gromacs/trajectoryanalysis/tests/refdata/SelectModuleTest_HandlesPDBOutputWithNonPDBInput.xml
@@ -0,0 +1,142 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <String Name="CommandLine">select -select 'resname RA RD and y &lt; 2.5'</String>
+  <OutputData Name="Data">
+    <AnalysisData Name="occupancy">
+      <DataFrame Name="Frame0">
+        <Real Name="X">1.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame1">
+        <Real Name="X">2.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame2">
+        <Real Name="X">3.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame3">
+        <Real Name="X">4.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame4">
+        <Real Name="X">5.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame5">
+        <Real Name="X">6.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame6">
+        <Real Name="X">7.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame7">
+        <Real Name="X">8.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame8">
+        <Real Name="X">9.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+    </AnalysisData>
+  </OutputData>
+  <OutputFiles Name="Files">
+    <String Name="-ofpdb"><![CDATA[
+TITLE     frame t= 0.000
+REMARK    THIS IS A SIMULATION BOX
+CRYST1  100.000  100.000  100.000  90.00  90.00  90.00 P 1           1
+MODEL        1
+ATOM      1  CB   RA     2      10.000  10.000   0.000  1.00  0.00            
+ATOM      2  S1   RA     2      10.000  20.000   0.000  0.50  0.00            
+ATOM      3  S2   RA     2      10.000  30.000   0.000  0.50  0.00            
+ATOM      4  CB   RB     3      10.000  40.000   0.000  0.00  0.00            
+ATOM      5  S1   RB     3      20.000  10.000   0.000  0.00  0.00            
+ATOM      6  S2   RB     3      20.000  20.000   0.000  0.00  0.00            
+ATOM      7  CB   RA     4      20.000  30.000   0.000  0.50  0.00            
+ATOM      8  S1   RA     4      20.000  40.000   0.000  0.50  0.00            
+ATOM      9  S2   RA     4      30.000  10.000   0.000  1.00  0.00            
+ATOM     10  CB   RC     5      30.000  20.000   0.000  0.00  0.00            
+ATOM     11  S1   RC     5      30.000  30.000   0.000  0.00  0.00            
+ATOM     12  S2   RC     5      30.000  40.000   0.000  0.00  0.00            
+ATOM     13  CB   RD     1      40.000  10.000   0.000  1.00  0.00            
+ATOM     14  S1   RD     1      40.000  20.000   0.000  0.50  0.00            
+ATOM     15  S2   RD     1      40.000  30.000   0.000  0.00  0.00            
+TER
+ENDMDL
+]]></String>
+  </OutputFiles>
+</ReferenceData>

diff --git a/src/gromacs/trajectoryanalysis/tests/refdata/SelectModuleTest_HandlesPDBOutputWithPDBInput.xml b/src/gromacs/trajectoryanalysis/tests/refdata/SelectModuleTest_HandlesPDBOutputWithPDBInput.xml
new file mode 100644 (file)
index 0000000..df81b35
--- /dev/null
+++ b/src/gromacs/trajectoryanalysis/tests/refdata/SelectModuleTest_HandlesPDBOutputWithPDBInput.xml
@@ -0,0 +1,142 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <String Name="CommandLine">select -select 'resname RA RD and y &lt; 2.5'</String>
+  <OutputData Name="Data">
+    <AnalysisData Name="occupancy">
+      <DataFrame Name="Frame0">
+        <Real Name="X">1.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame1">
+        <Real Name="X">2.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame2">
+        <Real Name="X">3.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame3">
+        <Real Name="X">4.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame4">
+        <Real Name="X">5.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame5">
+        <Real Name="X">6.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame6">
+        <Real Name="X">7.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame7">
+        <Real Name="X">8.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame8">
+        <Real Name="X">9.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+    </AnalysisData>
+  </OutputData>
+  <OutputFiles Name="Files">
+    <String Name="-ofpdb"><![CDATA[
+TITLE     frame t= 0.000
+REMARK    THIS IS A SIMULATION BOX
+CRYST1  100.000  100.000  100.000  90.00  90.00  90.00 P 1           1
+MODEL        1
+ATOM      1  CB   RA A   1      10.000  10.000   0.000  1.00  0.00            
+ATOM      2  S1   RA A   1      10.000  20.000   0.000  0.50  0.00            
+ATOM      3  S2   RA A   1      10.000  30.000   0.000  0.50  0.00            
+ATOM      4  CB   RB A   2      10.000  40.000   0.000  0.00  0.50            
+ATOM      5  S1   RB A   2      20.000  10.000   0.000  0.00  0.50            
+ATOM      6  S2   RB A   2      20.000  20.000   0.000  0.00  0.50            
+ATOM      7  CB   RA B   1      20.000  30.000   0.000  0.50  0.00            
+ATOM      8  S1   RA B   1      20.000  40.000   0.000  0.50  0.30            
+ATOM      9  S2   RA B   1      30.000  10.000   0.000  1.00  0.30            
+ATOM     10  CB   RC B   1A     30.000  20.000   0.000  0.00  0.00            
+ATOM     11  S1   RC B   1A     30.000  30.000   0.000  0.00  0.00            
+ATOM     12  S2   RC B   1A     30.000  40.000   0.000  0.00  0.30            
+ATOM     13  CB   RD B   2      40.000  10.000   0.000  1.00  0.00            
+ATOM     14  S1   RD B   2      40.000  20.000   0.000  0.50  0.00            
+ATOM     15  S2   RD B   2      40.000  30.000   0.000  0.00  0.00            
+TER
+ENDMDL
+]]></String>
+  </OutputFiles>
+</ReferenceData>

diff --git a/src/gromacs/trajectoryanalysis/tests/refdata/SelectModuleTest_HandlesSelectedPDBOutput.xml b/src/gromacs/trajectoryanalysis/tests/refdata/SelectModuleTest_HandlesSelectedPDBOutput.xml
new file mode 100644 (file)
index 0000000..0ef1feb
--- /dev/null
+++ b/src/gromacs/trajectoryanalysis/tests/refdata/SelectModuleTest_HandlesSelectedPDBOutput.xml
@@ -0,0 +1,135 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <String Name="CommandLine">select -select 'resname RA RD and y &lt; 2.5' -pdbatoms selected</String>
+  <OutputData Name="Data">
+    <AnalysisData Name="occupancy">
+      <DataFrame Name="Frame0">
+        <Real Name="X">1.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame1">
+        <Real Name="X">2.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame2">
+        <Real Name="X">3.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame3">
+        <Real Name="X">4.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame4">
+        <Real Name="X">5.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame5">
+        <Real Name="X">6.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame6">
+        <Real Name="X">7.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">1.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame7">
+        <Real Name="X">8.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.500000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+      <DataFrame Name="Frame8">
+        <Real Name="X">9.000000</Real>
+        <Sequence Name="Y">
+          <Int Name="Length">2</Int>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+          <DataValue>
+            <Real Name="Value">0.000000</Real>
+          </DataValue>
+        </Sequence>
+      </DataFrame>
+    </AnalysisData>
+  </OutputData>
+  <OutputFiles Name="Files">
+    <String Name="-ofpdb"><![CDATA[
+TITLE     frame t= 0.000
+REMARK    THIS IS A SIMULATION BOX
+CRYST1  100.000  100.000  100.000  90.00  90.00  90.00 P 1           1
+MODEL        1
+ATOM      1  CB   RA A   1      10.000  10.000   0.000  1.00  0.00            
+ATOM      2  S1   RA A   1      10.000  20.000   0.000  0.50  0.00            
+ATOM      3  S2   RA A   1      10.000  30.000   0.000  0.50  0.00            
+ATOM      7  CB   RA B   1      20.000  30.000   0.000  0.50  0.00            
+ATOM      8  S1   RA B   1      20.000  40.000   0.000  0.50  0.30            
+ATOM      9  S2   RA B   1      30.000  10.000   0.000  1.00  0.30            
+ATOM     13  CB   RD B   2      40.000  10.000   0.000  1.00  0.00            
+ATOM     14  S1   RD B   2      40.000  20.000   0.000  0.50  0.00            
+TER
+ENDMDL
+]]></String>
+  </OutputFiles>
+</ReferenceData>

diff --git a/src/gromacs/trajectoryanalysis/tests/select.cpp b/src/gromacs/trajectoryanalysis/tests/select.cpp
index 0d44ee835c57a03beb4136bb6097eb6d04796e1c..b358291035e86a097d4eec511aa06b77450fd4ad 100644 (file)
--- a/src/gromacs/trajectoryanalysis/tests/select.cpp
+++ b/src/gromacs/trajectoryanalysis/tests/select.cpp
@@ -70,12 +70,67 @@ TEST_F(SelectModuleTest, BasicTest)
         "-select", "y < 2.5", "resname RA"
     };
     setTopology("simple.gro");
+    setTrajectory("simple.gro");
     setOutputFile("-oi", "index.dat");
     setOutputFile("-on", "index.ndx");
     excludeDataset("cfrac");
     runTest(CommandLine::create(cmdline));
 }
 
+TEST_F(SelectModuleTest, HandlesPDBOutputWithNonPDBInput)
+{
+    const char *const cmdline[] = {
+        "select",
+        "-select", "resname RA RD and y < 2.5"
+    };
+    setTopology("simple.gro");
+    setTrajectory("simple.gro");
+    includeDataset("occupancy");
+    setOutputFile("-ofpdb", "occupancy.pdb");
+    runTest(CommandLine::create(cmdline));
+}
+
+TEST_F(SelectModuleTest, HandlesPDBOutputWithPDBInput)
+{
+    const char *const cmdline[] = {
+        "select",
+        "-select", "resname RA RD and y < 2.5"
+    };
+    setTopology("simple.pdb");
+    setTrajectory("simple.gro");
+    includeDataset("occupancy");
+    setOutputFile("-ofpdb", "occupancy.pdb");
+    runTest(CommandLine::create(cmdline));
+}
+
+TEST_F(SelectModuleTest, HandlesMaxPDBOutput)
+{
+    const char *const cmdline[] = {
+        "select",
+        "-select", "resname RA RD and y < 2.5",
+        "-pdbatoms", "maxsel"
+    };
+    setTopology("simple.pdb");
+    setTrajectory("simple.gro");
+    includeDataset("occupancy");
+    setOutputFile("-ofpdb", "occupancy.pdb");
+    runTest(CommandLine::create(cmdline));
+}
+
+TEST_F(SelectModuleTest, HandlesSelectedPDBOutput)
+{
+    const char *const cmdline[] = {
+        "select",
+        "-select", "resname RA RD and y < 2.5",
+        "-pdbatoms", "selected"
+    };
+    setTopology("simple.pdb");
+    setTrajectory("simple.gro");
+    includeDataset("occupancy");
+    setOutputFile("-ofpdb", "occupancy.pdb");
+    runTest(CommandLine::create(cmdline));
+}
+
 TEST_F(SelectModuleTest, HandlesDumpOption)
 {
     const char *const cmdline[] = {

diff --git a/src/gromacs/trajectoryanalysis/tests/simple.gro b/src/gromacs/trajectoryanalysis/tests/simple.gro
index 13e93cd2986cdbbfa327f330acab5674df564b02..d8a6579e1dfc520f875bc37b0605c4248e3aac96 100644 (file)
--- a/src/gromacs/trajectoryanalysis/tests/simple.gro
+++ b/src/gromacs/trajectoryanalysis/tests/simple.gro
@@ -16,3 +16,21 @@ Test system
     1RD      S1   14   4.000   2.000   0.000
     1RD      S2   15   4.000   3.000   0.000
   10.00000  10.00000  10.00000
+Test system
+ 15
+    2RA      CB    1   1.000   2.000   0.000
+    2RA      S1    2   1.000   3.000   0.000
+    2RA      S2    3   1.000   1.000   0.000
+    3RB      CB    4   1.000   2.000   0.000
+    3RB      S1    5   2.000   3.000   0.000
+    3RB      S2    6   2.000   4.000   0.000
+    4RA      CB    7   2.000   1.000   0.000
+    4RA      S1    8   2.000   2.000   0.000
+    4RA      S2    9   3.000   2.000   0.000
+    5RC      CB   10   3.000   3.000   0.000
+    5RC      S1   11   3.000   4.000   0.000
+    5RC      S2   12   3.000   1.000   0.000
+    1RD      CB   13   4.000   2.000   0.000
+    1RD      S1   14   4.000   3.000   0.000
+    1RD      S2   15   4.000   4.000   0.000
+  10.00000  10.00000  10.00000

diff --git a/src/gromacs/trajectoryanalysis/tests/simple.pdb b/src/gromacs/trajectoryanalysis/tests/simple.pdb
new file mode 100644 (file)
index 0000000..156b948
--- /dev/null
+++ b/src/gromacs/trajectoryanalysis/tests/simple.pdb
@@ -0,0 +1,20 @@
+TITLE     Test system t=   0.00000
+CRYST1  100.000  100.000  100.000  90.00  90.00  90.00 P 1           1
+MODEL        1
+ATOM      1  CB   RA A   1      10.000  10.000   0.000  1.00  0.00            
+ATOM      2  S1   RA A   1      10.000  20.000   0.000  1.00  0.00            
+ATOM      3  S2   RA A   1      10.000  30.000   0.000  1.00  0.00            
+ATOM      4  CB   RB A   2      10.000  40.000   0.000  1.00  0.50            
+ATOM      5  S1   RB A   2      20.000  10.000   0.000  1.00  0.50            
+ATOM      6  S2   RB A   2      20.000  20.000   0.000  1.00  0.50            
+ATOM      7  CB   RA B   1      20.000  30.000   0.000  1.00  0.00            
+ATOM      8  S1   RA B   1      20.000  40.000   0.000  1.00  0.30            
+ATOM      9  S2   RA B   1      30.000  10.000   0.000  1.00  0.30            
+ATOM     10  CB   RC B   1A     30.000  20.000   0.000  1.00  0.00            
+ATOM     11  S1   RC B   1A     30.000  30.000   0.000  1.00  0.00            
+ATOM     12  S2   RC B   1A     30.000  40.000   0.000  1.00  0.30            
+ATOM     13  CB A RD B   2      40.000  10.000   0.000  1.00  0.00            
+ATOM     14  S1 A RD B   2      40.000  20.000   0.000  1.00  0.00            
+ATOM     15  S2 A RD B   2      40.000  30.000   0.000  1.00  0.00            
+TER
+ENDMDL

diff --git a/src/gromacs/utility/arrayref.h b/src/gromacs/utility/arrayref.h
index e38fa5987df6b864cd87ae4bd61a54f40af97d53..471839de68ff68fedf4a85bc9ef95b67289b3001 100644 (file)
--- a/src/gromacs/utility/arrayref.h
+++ b/src/gromacs/utility/arrayref.h
@@ -175,6 +175,9 @@ class ConstArrayRef
         //! Returns the last element in the container.
         const_reference back() const { return *(end_ - 1); }
 
+        //! Returns a raw pointer to the contents of the array.
+        const_pointer data() const { return begin_; }
+
         /*! \brief
          * Swaps referenced memory with the other object.
          *