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, 6 months 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
#	User	Rev	Content
1	gezelter	1865	/*
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	gezelter	1876	#ifdef _MSC_VER
47			#define isnan(x) _isnan((x))
48			#define isinf(x) (!_finite(x) && !_isnan(x))
49			#endif
50	gezelter	1865
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	gezelter	1881	}
145	gezelter	1865	}
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	gezelter	1979	outStream << "# 95% confidence intervals in those quantities follow:\n";
191	gezelter	1865	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	gezelter	1979	writeErrorBars( outStream, outputData, j );
204	gezelter	1865	}
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	gezelter	1875	RealType r;
230	gezelter	1865	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	gezelter	1875	"\t%s for bin %u",
235	gezelter	1865	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	gezelter	1875	Vector3d v;
244	gezelter	1865	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	gezelter	1875	"\t%s for bin %u",
251	gezelter	1865	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	gezelter	1979	void SpatialStatistics::writeErrorBars(ostream& os, OutputData* dat,
261	gezelter	1865	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	gezelter	1875	RealType r;
268	gezelter	1979	dynamic_cast<Accumulator*>(dat->accumulator[bin])->get95percentConfidenceInterval(r);
269	gezelter	1865	if (isinf(r) \|\| isnan(r) ) {
270			sprintf( painCave.errMsg,
271			"SpatialStatistics detected a numerical error writing:\n"
272	gezelter	1875	"\tstandard deviation of %s for bin %u",
273	gezelter	1865	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	gezelter	1875	Vector3d v;
282	gezelter	1979	dynamic_cast<VectorAccumulator*>(dat->accumulator[bin])->get95percentConfidenceInterval(v);
283	gezelter	1865	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	gezelter	1875	"\tstandard deviation of %s for bin %u",
289	gezelter	1865	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	gezelter	1874	SlabStatistics::~SlabStatistics() {
326			}
327
328
329	gezelter	1865	void SlabStatistics::processFrame(int istep) {
330			RealType z;
331	gezelter	1887
332	gezelter	1865	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	gezelter	1875	SpatialStatistics(info, filename, sele, nbins), coordinateOrigin_(V3Zero) {
355	gezelter	1865
356			binWidth_ = 1.0;
357	gezelter	1945
358			Globals* simParams = info->getSimParams();
359			RNEMDParameters* rnemdParams = simParams->getRNEMDParameters();
360			bool hasCoordinateOrigin = rnemdParams->haveCoordinateOrigin();
361	gezelter	1875
362	gezelter	1945	if (hasCoordinateOrigin) {
363			coordinateOrigin_ = rnemdParams->getCoordinateOrigin();
364			} else {
365			coordinateOrigin_ = V3Zero;
366			}
367
368	gezelter	1865	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	gezelter	1874	ShellStatistics::~ShellStatistics() {
385			}
386
387	gezelter	1945	int ShellStatistics::getBin(Vector3d pos) {
388	gezelter	1865	Vector3d rPos = pos - coordinateOrigin_;
389			return int(rPos.length() / binWidth_);
390			}
391			}
392