applications/staticProps/SpatialStatistics.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
 * [4] Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). 
 */

#include "applications/staticProps/SpatialStatistics.hpp"
#include "io/DumpReader.hpp"
#include "primitives/Molecule.hpp"
#ifdef _MSC_VER
#define isnan(x) _isnan((x))
#define isinf(x) (!_finite(x) && !_isnan(x))
#endif

namespace OpenMD {
  
  SpatialStatistics::SpatialStatistics(SimInfo* info, const string& filename, 
                                       const string& sele, int nbins)
    : StaticAnalyser(info, filename), selectionScript_(sele),  evaluator_(info),
      seleMan_(info), nBins_(nbins){
    
    evaluator_.loadScriptString(sele);
    if (!evaluator_.isDynamic()) {
      seleMan_.setSelectionSet(evaluator_.evaluate());
    }
    
    // Pre-load an OutputData for the count of objects:
    counts_ = new OutputData;
    counts_->units =  "objects";
    counts_->title =  "Objects";
    counts_->dataType = odtReal;
    counts_->dataHandling = odhTotal;
    counts_->accumulator.reserve(nBins_);
    for (int i = 0; i < nBins_; i++) 
      counts_->accumulator.push_back( new Accumulator() );
    
    setOutputName(getPrefix(filename) + ".spst");
  }

  SpatialStatistics::~SpatialStatistics() {
    vector<OutputData*>::iterator i;
    OutputData* outputData;
    
    for(outputData = beginOutputData(i); outputData; 
        outputData = nextOutputData(i)) {
      delete outputData;
    }
    data_.clear();

    delete counts_;
  }


  void SpatialStatistics::process() {

    DumpReader reader(info_, dumpFilename_);    
    int nFrames = reader.getNFrames();
    nProcessed_ = nFrames/step_;

    for (int istep = 0; istep < nFrames; istep += step_) {
      reader.readFrame(istep);
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
      processFrame(istep);
    }
    writeOutput();
  }


  void SpatialStatistics::processFrame(int istep) {
    Molecule* mol;
    RigidBody* rb;
    StuntDouble* sd;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
    int i;
    
    for (mol = info_->beginMolecule(mi); mol != NULL; 
         mol = info_->nextMolecule(mi)) {
      
      // change the positions of atoms which belong to the rigidbodies
      
      for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
           rb = mol->nextRigidBody(rbIter)) {
        rb->updateAtoms();
      }
    }
    
    if (evaluator_.isDynamic()) {
      seleMan_.setSelectionSet(evaluator_.evaluate());
    }
    
    // loop over the selected atoms:
    
    for (sd = seleMan_.beginSelected(i); sd != NULL; 
         sd = seleMan_.nextSelected(i)) {
      
      // figure out where that object is:
      
      Vector3d pos = sd->getPos();
      
      int bin = getBin(pos);
      
      // forward the work of statistics on to the subclass:
      
      processStuntDouble( sd, bin );

      dynamic_cast<Accumulator *>(counts_->accumulator[bin])->add(1);
    }
  }
  

  void SpatialStatistics::writeOutput() {
    
    vector<OutputData*>::iterator i;
    OutputData* outputData;
    
    ofstream outStream(outputFilename_.c_str());
    if (outStream.is_open()) {
      
      //write title
      outStream << "# SPATIAL STATISTICS\n";
      outStream << "#";
      
      for(outputData = beginOutputData(i); outputData; 
          outputData = nextOutputData(i)) {
        outStream << "\t" << outputData->title << 
          "(" << outputData->units << ")";
        // add some extra tabs for column alignment
        if (outputData->dataType == odtVector3) outStream << "\t\t";
      }
      
      outStream << std::endl;
      
      outStream.precision(8);
      
      for (int j = 0; j < nBins_; j++) {        
        
        int counts = counts_->accumulator[j]->count();

        if (counts > 0) {
          for(outputData = beginOutputData(i); outputData; 
              outputData = nextOutputData(i)) {
            
            int n = outputData->accumulator[j]->count();
            if (n != 0) {
              writeData( outStream, outputData, j );
            }
          }
          outStream << std::endl;
        }
      }
        
      outStream << "#######################################################\n";
      outStream << "# 95% confidence intervals in those quantities follow:\n";
      outStream << "#######################################################\n";
      
      for (int j = 0; j < nBins_; j++) {
        int counts = counts_->accumulator[j]->count();
        if (counts > 0) {
          
          outStream << "#";
          for(outputData = beginOutputData(i); outputData; 
              outputData = nextOutputData(i)) {
            
            int n = outputData->accumulator[j]->count();
            if (n != 0) {
              writeErrorBars( outStream, outputData, j );
            }
          }
          outStream << std::endl;
        }
      }
      
      outStream.flush();
      outStream.close();      
      
    } else {      
      sprintf(painCave.errMsg, "SpatialStatistics: unable to open %s\n", 
              outputFilename_.c_str());
      painCave.isFatal = 1;
      simError();  
    }   
  }
  
  
  void SpatialStatistics::writeData(ostream& os, OutputData* dat, 
                                    unsigned int bin) {
    assert(int(bin) < nBins_);
    int n = dat->accumulator[bin]->count();
    if (n == 0) return;

    if( dat->dataType == odtReal ) {
      RealType r;
      dynamic_cast<Accumulator*>(dat->accumulator[bin])->getAverage(r);      
      if (isinf(r) || isnan(r) ) {      
        sprintf( painCave.errMsg,
                 "SpatialStatistics detected a numerical error writing:\n"
                 "\t%s for bin %u",
                 dat->title.c_str(), bin);
        painCave.isFatal = 1;
        simError();
      }
      if (dat->dataHandling == odhTotal) r *= dat->accumulator[bin]->count();
      os << "\t" << r;      

    } else if ( dat->dataType == odtVector3 ) {
      Vector3d v;
      dynamic_cast<VectorAccumulator*>(dat->accumulator[bin])->getAverage(v);
      if (isinf(v[0]) || isnan(v[0]) || 
          isinf(v[1]) || isnan(v[1]) || 
          isinf(v[2]) || isnan(v[2]) ) {      
        sprintf( painCave.errMsg,
                 "SpatialStatistics detected a numerical error writing:\n"
                 "\t%s for bin %u",
                 dat->title.c_str(), bin);
        painCave.isFatal = 1;
        simError();
      }
      if (dat->dataHandling == odhTotal) v *= dat->accumulator[bin]->count();
      os << "\t" << v[0] << "\t" << v[1] << "\t" << v[2];
    }
  }

  void SpatialStatistics::writeErrorBars(ostream& os, OutputData* dat, 
                                    unsigned int bin) {
    assert(int(bin) < nBins_);
    int n = dat->accumulator[bin]->count();
    if (n == 0) return;

    if( dat->dataType == odtReal ) {
      RealType r;
      dynamic_cast<Accumulator*>(dat->accumulator[bin])->get95percentConfidenceInterval(r);      
      if (isinf(r) || isnan(r) ) {      
        sprintf( painCave.errMsg,
                 "SpatialStatistics detected a numerical error writing:\n"
                 "\tstandard deviation of %s for bin %u",
                 dat->title.c_str(), bin);
        painCave.isFatal = 1;
        simError();
      }
      if (dat->dataHandling == odhTotal) r *= dat->accumulator[bin]->count();
      os << "\t" << r;      

    } else if ( dat->dataType == odtVector3 ) {
      Vector3d v;
      dynamic_cast<VectorAccumulator*>(dat->accumulator[bin])->get95percentConfidenceInterval(v);
      if (isinf(v[0]) || isnan(v[0]) || 
          isinf(v[1]) || isnan(v[1]) || 
          isinf(v[2]) || isnan(v[2]) ) {      
        sprintf( painCave.errMsg,
                 "SpatialStatistics detected a numerical error writing:\n"
                 "\tstandard deviation of %s for bin %u",
                 dat->title.c_str(), bin);
        painCave.isFatal = 1;
        simError();
      }
      if (dat->dataHandling == odhTotal) v *= dat->accumulator[bin]->count();
      os << "\t" << v[0] << "\t" << v[1] << "\t" << v[2];
    }
  }
  
  
  OutputData* SpatialStatistics::beginOutputData(vector<OutputData*>::iterator& i) {
    i = data_.begin();
    return i != data_.end()? *i : NULL;
  }

  OutputData* SpatialStatistics::nextOutputData(vector<OutputData*>::iterator& i){
    ++i;
    return i != data_.end()? *i: NULL;
  }


  SlabStatistics::SlabStatistics(SimInfo* info, const string& filename, 
                                 const string& sele, int nbins) : 
    SpatialStatistics(info, filename, sele, nbins) {
    
    z_ = new OutputData;
    z_->units =  "Angstroms";
    z_->title =  "Z";
    z_->dataType = odtReal;
    z_->dataHandling = odhAverage;
    z_->accumulator.reserve(nbins);
    for (int i = 0; i < nbins; i++) 
      z_->accumulator.push_back( new Accumulator() );
    data_.push_back(z_);
  }

  SlabStatistics::~SlabStatistics() {
  }


  void SlabStatistics::processFrame(int istep) {
    RealType z;

    hmat_ = currentSnapshot_->getHmat();
    for (int i = 0; i < nBins_; i++) {
      z = (((RealType)i + 0.5) / (RealType)nBins_) * hmat_(2,2);
      dynamic_cast<Accumulator*>(z_->accumulator[i])->add(z);
    }
    volume_ = currentSnapshot_->getVolume();

    SpatialStatistics::processFrame(istep);
  }

  int SlabStatistics::getBin(Vector3d pos) {
    currentSnapshot_->wrapVector(pos);
    // which bin is this stuntdouble in?
    // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
    // Shift molecules by half a box to have bins start at 0
    // The modulo operator is used to wrap the case when we are 
    // beyond the end of the bins back to the beginning.
    return int(nBins_ * (pos.z() / hmat_(2,2) + 0.5)) % nBins_;  
  }

  ShellStatistics::ShellStatistics(SimInfo* info, const string& filename, 
                                   const string& sele, int nbins) : 
    SpatialStatistics(info, filename, sele, nbins), coordinateOrigin_(V3Zero) {
    
    binWidth_ = 1.0;

    Globals* simParams = info->getSimParams();
    RNEMDParameters* rnemdParams = simParams->getRNEMDParameters();
    bool hasCoordinateOrigin = rnemdParams->haveCoordinateOrigin();
    
    if (hasCoordinateOrigin) {
      coordinateOrigin_ = rnemdParams->getCoordinateOrigin();
    } else {
      coordinateOrigin_ = V3Zero;
    }
    
    r_ = new OutputData;
    r_->units =  "Angstroms";
    r_->title =  "R";
    r_->dataType = odtReal;
    r_->dataHandling = odhAverage;
    r_->accumulator.reserve(nbins);
    for (int i = 0; i < nbins; i++)
      r_->accumulator.push_back( new Accumulator() );    
    data_.push_back(r_);

    for (int i = 0; i < nbins; i++) {
      RealType r = (((RealType)i + 0.5) * binWidth_);
      dynamic_cast<Accumulator*>(r_->accumulator[i])->add(r);
    }
  }

  ShellStatistics::~ShellStatistics() {
  }

  int ShellStatistics::getBin(Vector3d pos) { 
    Vector3d rPos = pos - coordinateOrigin_;
    return int(rPos.length() / binWidth_);
  }
}

Revision:	1979
Committed:	Sat Apr 5 20:56:01 2014 UTC (11 years ago) by gezelter
File size:	13049 byte(s)
Log Message:	* Changed the stdDev printouts in RNEMD stats to the 95% confidence intervals * Added the ability to specify non-bonded constraints in a molecule * Added the ability to do selection offsets in the pAngle staticProps module
#	Content
1	/*
2	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	*
4	* The University of Notre Dame grants you ("Licensee") a
5	* non-exclusive, royalty free, license to use, modify and
6	* redistribute this software in source and binary code form, provided
7	* that the following conditions are met:
8	*
9	* 1. Redistributions of source code must retain the above copyright
10	* notice, this list of conditions and the following disclaimer.
11	*
12	* 2. Redistributions in binary form must reproduce the above copyright
13	* notice, this list of conditions and the following disclaimer in the
14	* documentation and/or other materials provided with the
15	* distribution.
16	*
17	* This software is provided "AS IS," without a warranty of any
18	* kind. All express or implied conditions, representations and
19	* warranties, including any implied warranty of merchantability,
20	* fitness for a particular purpose or non-infringement, are hereby
21	* excluded. The University of Notre Dame and its licensors shall not
22	* be liable for any damages suffered by licensee as a result of
23	* using, modifying or distributing the software or its
24	* derivatives. In no event will the University of Notre Dame or its
25	* licensors be liable for any lost revenue, profit or data, or for
26	* direct, indirect, special, consequential, incidental or punitive
27	* damages, however caused and regardless of the theory of liability,
28	* arising out of the use of or inability to use software, even if the
29	* University of Notre Dame has been advised of the possibility of
30	* such damages.
31	*
32	* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33	* research, please cite the appropriate papers when you publish your
34	* work. Good starting points are:
35	*
36	* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37	* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40	* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	*/
42
43	#include "applications/staticProps/SpatialStatistics.hpp"
44	#include "io/DumpReader.hpp"
45	#include "primitives/Molecule.hpp"
46	#ifdef _MSC_VER
47	#define isnan(x) _isnan((x))
48	#define isinf(x) (!_finite(x) && !_isnan(x))
49	#endif
50
51	namespace OpenMD {
52
53	SpatialStatistics::SpatialStatistics(SimInfo* info, const string& filename,
54	const string& sele, int nbins)
55	: StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info),
56	seleMan_(info), nBins_(nbins){
57
58	evaluator_.loadScriptString(sele);
59	if (!evaluator_.isDynamic()) {
60	seleMan_.setSelectionSet(evaluator_.evaluate());
61	}
62
63	// Pre-load an OutputData for the count of objects:
64	counts_ = new OutputData;
65	counts_->units = "objects";
66	counts_->title = "Objects";
67	counts_->dataType = odtReal;
68	counts_->dataHandling = odhTotal;
69	counts_->accumulator.reserve(nBins_);
70	for (int i = 0; i < nBins_; i++)
71	counts_->accumulator.push_back( new Accumulator() );
72
73	setOutputName(getPrefix(filename) + ".spst");
74	}
75
76	SpatialStatistics::~SpatialStatistics() {
77	vector<OutputData*>::iterator i;
78	OutputData* outputData;
79
80	for(outputData = beginOutputData(i); outputData;
81	outputData = nextOutputData(i)) {
82	delete outputData;
83	}
84	data_.clear();
85
86	delete counts_;
87	}
88
89
90	void SpatialStatistics::process() {
91
92	DumpReader reader(info_, dumpFilename_);
93	int nFrames = reader.getNFrames();
94	nProcessed_ = nFrames/step_;
95
96	for (int istep = 0; istep < nFrames; istep += step_) {
97	reader.readFrame(istep);
98	currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
99	processFrame(istep);
100	}
101	writeOutput();
102	}
103
104
105	void SpatialStatistics::processFrame(int istep) {
106	Molecule* mol;
107	RigidBody* rb;
108	StuntDouble* sd;
109	SimInfo::MoleculeIterator mi;
110	Molecule::RigidBodyIterator rbIter;
111	int i;
112
113	for (mol = info_->beginMolecule(mi); mol != NULL;
114	mol = info_->nextMolecule(mi)) {
115
116	// change the positions of atoms which belong to the rigidbodies
117
118	for (rb = mol->beginRigidBody(rbIter); rb != NULL;
119	rb = mol->nextRigidBody(rbIter)) {
120	rb->updateAtoms();
121	}
122	}
123
124	if (evaluator_.isDynamic()) {
125	seleMan_.setSelectionSet(evaluator_.evaluate());
126	}
127
128	// loop over the selected atoms:
129
130	for (sd = seleMan_.beginSelected(i); sd != NULL;
131	sd = seleMan_.nextSelected(i)) {
132
133	// figure out where that object is:
134
135	Vector3d pos = sd->getPos();
136
137	int bin = getBin(pos);
138
139	// forward the work of statistics on to the subclass:
140
141	processStuntDouble( sd, bin );
142
143	dynamic_cast<Accumulator *>(counts_->accumulator[bin])->add(1);
144	}
145	}
146
147
148	void SpatialStatistics::writeOutput() {
149
150	vector<OutputData*>::iterator i;
151	OutputData* outputData;
152
153	ofstream outStream(outputFilename_.c_str());
154	if (outStream.is_open()) {
155
156	//write title
157	outStream << "# SPATIAL STATISTICS\n";
158	outStream << "#";
159
160	for(outputData = beginOutputData(i); outputData;
161	outputData = nextOutputData(i)) {
162	outStream << "\t" << outputData->title <<
163	"(" << outputData->units << ")";
164	// add some extra tabs for column alignment
165	if (outputData->dataType == odtVector3) outStream << "\t\t";
166	}
167
168	outStream << std::endl;
169
170	outStream.precision(8);
171
172	for (int j = 0; j < nBins_; j++) {
173
174	int counts = counts_->accumulator[j]->count();
175
176	if (counts > 0) {
177	for(outputData = beginOutputData(i); outputData;
178	outputData = nextOutputData(i)) {
179
180	int n = outputData->accumulator[j]->count();
181	if (n != 0) {
182	writeData( outStream, outputData, j );
183	}
184	}
185	outStream << std::endl;
186	}
187	}
188
189	outStream << "#######################################################\n";
190	outStream << "# 95% confidence intervals in those quantities follow:\n";
191	outStream << "#######################################################\n";
192
193	for (int j = 0; j < nBins_; j++) {
194	int counts = counts_->accumulator[j]->count();
195	if (counts > 0) {
196
197	outStream << "#";
198	for(outputData = beginOutputData(i); outputData;
199	outputData = nextOutputData(i)) {
200
201	int n = outputData->accumulator[j]->count();
202	if (n != 0) {
203	writeErrorBars( outStream, outputData, j );
204	}
205	}
206	outStream << std::endl;
207	}
208	}
209
210	outStream.flush();
211	outStream.close();
212
213	} else {
214	sprintf(painCave.errMsg, "SpatialStatistics: unable to open %s\n",
215	outputFilename_.c_str());
216	painCave.isFatal = 1;
217	simError();
218	}
219	}
220
221
222	void SpatialStatistics::writeData(ostream& os, OutputData* dat,
223	unsigned int bin) {
224	assert(int(bin) < nBins_);
225	int n = dat->accumulator[bin]->count();
226	if (n == 0) return;
227
228	if( dat->dataType == odtReal ) {
229	RealType r;
230	dynamic_cast<Accumulator*>(dat->accumulator[bin])->getAverage(r);
231	if (isinf(r) \|\| isnan(r) ) {
232	sprintf( painCave.errMsg,
233	"SpatialStatistics detected a numerical error writing:\n"
234	"\t%s for bin %u",
235	dat->title.c_str(), bin);
236	painCave.isFatal = 1;
237	simError();
238	}
239	if (dat->dataHandling == odhTotal) r *= dat->accumulator[bin]->count();
240	os << "\t" << r;
241
242	} else if ( dat->dataType == odtVector3 ) {
243	Vector3d v;
244	dynamic_cast<VectorAccumulator*>(dat->accumulator[bin])->getAverage(v);
245	if (isinf(v[0]) \|\| isnan(v[0]) \|\|
246	isinf(v[1]) \|\| isnan(v[1]) \|\|
247	isinf(v[2]) \|\| isnan(v[2]) ) {
248	sprintf( painCave.errMsg,
249	"SpatialStatistics detected a numerical error writing:\n"
250	"\t%s for bin %u",
251	dat->title.c_str(), bin);
252	painCave.isFatal = 1;
253	simError();
254	}
255	if (dat->dataHandling == odhTotal) v *= dat->accumulator[bin]->count();
256	os << "\t" << v[0] << "\t" << v[1] << "\t" << v[2];
257	}
258	}
259
260	void SpatialStatistics::writeErrorBars(ostream& os, OutputData* dat,
261	unsigned int bin) {
262	assert(int(bin) < nBins_);
263	int n = dat->accumulator[bin]->count();
264	if (n == 0) return;
265
266	if( dat->dataType == odtReal ) {
267	RealType r;
268	dynamic_cast<Accumulator*>(dat->accumulator[bin])->get95percentConfidenceInterval(r);
269	if (isinf(r) \|\| isnan(r) ) {
270	sprintf( painCave.errMsg,
271	"SpatialStatistics detected a numerical error writing:\n"
272	"\tstandard deviation of %s for bin %u",
273	dat->title.c_str(), bin);
274	painCave.isFatal = 1;
275	simError();
276	}
277	if (dat->dataHandling == odhTotal) r *= dat->accumulator[bin]->count();
278	os << "\t" << r;
279
280	} else if ( dat->dataType == odtVector3 ) {
281	Vector3d v;
282	dynamic_cast<VectorAccumulator*>(dat->accumulator[bin])->get95percentConfidenceInterval(v);
283	if (isinf(v[0]) \|\| isnan(v[0]) \|\|
284	isinf(v[1]) \|\| isnan(v[1]) \|\|
285	isinf(v[2]) \|\| isnan(v[2]) ) {
286	sprintf( painCave.errMsg,
287	"SpatialStatistics detected a numerical error writing:\n"
288	"\tstandard deviation of %s for bin %u",
289	dat->title.c_str(), bin);
290	painCave.isFatal = 1;
291	simError();
292	}
293	if (dat->dataHandling == odhTotal) v *= dat->accumulator[bin]->count();
294	os << "\t" << v[0] << "\t" << v[1] << "\t" << v[2];
295	}
296	}
297
298
299	OutputData* SpatialStatistics::beginOutputData(vector<OutputData*>::iterator& i) {
300	i = data_.begin();
301	return i != data_.end()? *i : NULL;
302	}
303
304	OutputData* SpatialStatistics::nextOutputData(vector<OutputData*>::iterator& i){
305	++i;
306	return i != data_.end()? *i: NULL;
307	}
308
309
310	SlabStatistics::SlabStatistics(SimInfo* info, const string& filename,
311	const string& sele, int nbins) :
312	SpatialStatistics(info, filename, sele, nbins) {
313
314	z_ = new OutputData;
315	z_->units = "Angstroms";
316	z_->title = "Z";
317	z_->dataType = odtReal;
318	z_->dataHandling = odhAverage;
319	z_->accumulator.reserve(nbins);
320	for (int i = 0; i < nbins; i++)
321	z_->accumulator.push_back( new Accumulator() );
322	data_.push_back(z_);
323	}
324
325	SlabStatistics::~SlabStatistics() {
326	}
327
328
329	void SlabStatistics::processFrame(int istep) {
330	RealType z;
331
332	hmat_ = currentSnapshot_->getHmat();
333	for (int i = 0; i < nBins_; i++) {
334	z = (((RealType)i + 0.5) / (RealType)nBins_) * hmat_(2,2);
335	dynamic_cast<Accumulator*>(z_->accumulator[i])->add(z);
336	}
337	volume_ = currentSnapshot_->getVolume();
338
339	SpatialStatistics::processFrame(istep);
340	}
341
342	int SlabStatistics::getBin(Vector3d pos) {
343	currentSnapshot_->wrapVector(pos);
344	// which bin is this stuntdouble in?
345	// wrapped positions are in the range [-0.5hmat(2,2), +0.5hmat(2,2)]
346	// Shift molecules by half a box to have bins start at 0
347	// The modulo operator is used to wrap the case when we are
348	// beyond the end of the bins back to the beginning.
349	return int(nBins_ * (pos.z() / hmat_(2,2) + 0.5)) % nBins_;
350	}
351
352	ShellStatistics::ShellStatistics(SimInfo* info, const string& filename,
353	const string& sele, int nbins) :
354	SpatialStatistics(info, filename, sele, nbins), coordinateOrigin_(V3Zero) {
355
356	binWidth_ = 1.0;
357
358	Globals* simParams = info->getSimParams();
359	RNEMDParameters* rnemdParams = simParams->getRNEMDParameters();
360	bool hasCoordinateOrigin = rnemdParams->haveCoordinateOrigin();
361
362	if (hasCoordinateOrigin) {
363	coordinateOrigin_ = rnemdParams->getCoordinateOrigin();
364	} else {
365	coordinateOrigin_ = V3Zero;
366	}
367
368	r_ = new OutputData;
369	r_->units = "Angstroms";
370	r_->title = "R";
371	r_->dataType = odtReal;
372	r_->dataHandling = odhAverage;
373	r_->accumulator.reserve(nbins);
374	for (int i = 0; i < nbins; i++)
375	r_->accumulator.push_back( new Accumulator() );
376	data_.push_back(r_);
377
378	for (int i = 0; i < nbins; i++) {
379	RealType r = (((RealType)i + 0.5) * binWidth_);
380	dynamic_cast<Accumulator*>(r_->accumulator[i])->add(r);
381	}
382	}
383
384	ShellStatistics::~ShellStatistics() {
385	}
386
387	int ShellStatistics::getBin(Vector3d pos) {
388	Vector3d rPos = pos - coordinateOrigin_;
389	return int(rPos.length() / binWidth_);
390	}
391	}
392