applications/staticProps/TetrahedralityParamZ.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, 24107 (2008).          
 * [4] Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
 *  Created by J. Daniel Gezelter on 09/26/06.
 *  @author  J. Daniel Gezelter
 *  @version $Id: BondOrderParameter.cpp 1442 2010-05-10 17:28:26Z gezelter $
 *
 */
 
#include "applications/staticProps/TetrahedralityParamZ.hpp"
#include "utils/simError.h"
#include "io/DumpReader.hpp"
#include "primitives/Molecule.hpp"
#include "utils/NumericConstant.hpp"
#include <vector>
#include <algorithm>
#include <fstream>

using namespace std;

namespace OpenMD 
{
  TetrahedralityParamZ::TetrahedralityParamZ(SimInfo* info, 
                                           const std::string& filename, 
                                           const std::string& sele,
                                           double rCut, int nzbins) : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan1_(info),seleMan2_(info), nZBins_(nzbins)
  {
    //nZBins_ = 50;
    //std ::cerrnZBins_:"<<nZBins_<<"\t"<<"nzbins:"<<nzbins<<endl;
    // nZBins_ = 90;
    //fixed numbe of bins
    count_.resize(nZBins_);
    sliceSDLists_.resize(nZBins_);
    Qave_.resize(nZBins_);

    setOutputName(getPrefix(filename) + ".q");
    
    evaluator_.loadScriptString(sele);
    if (!evaluator_.isDynamic()) 
      {
        seleMan1_.setSelectionSet(evaluator_.evaluate());
        seleMan2_.setSelectionSet(evaluator_.evaluate());
      }

    // Set up cutoff radius:
    rCut_ = rCut;

    // Q can take values from 0 to 1
    MinQ_ = 0.0;
    MaxQ_ = 1.1;
    deltaQ_ = (MaxQ_ - MinQ_)/nzbins;
  }

  TetrahedralityParamZ::~TetrahedralityParamZ() 
  {
    Q_histogram_.clear(); 
  }
  
  void TetrahedralityParamZ::initalizeHistogram() 
  {
    std::fill(Q_histogram_.begin(), Q_histogram_.end(), 0);
  }
  

  void TetrahedralityParamZ::process() 
  {
    Molecule* mol;
    StuntDouble* sd;
    StuntDouble* sd2;
    StuntDouble* sdi;
    StuntDouble* sdj;
    RigidBody* rb;
    int myIndex;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
    Vector3d vec;
    Vector3d ri, rj, rk, rik, rkj, dposition, tposition;
    RealType r;
    RealType cospsi;
    RealType Qk;

    std::vector<std::pair<RealType,StuntDouble*> > myNeighbors;
    int isd1, isd2;
    cerr << "After Creation of variables in TP:process()\n";
    DumpReader reader(info_, dumpFilename_);    
    cerr << "The DumpReader was created?\n";
    cerr << "nZbins: " << nZBins_ << "\n";
    int nFrames = reader.getNFrames();
    frameCounter_ = 0;
    nProcessed_=nFrames/step_;
    reader.readFrame(0); 
    currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
    Mat3x3d hmat = currentSnapshot_->getHmat();
    zBox_.push_back(hmat(2,2));
    
    RealType halfBoxZ_ = hmat(2,2) / 2.0; 

    Distorted_.clear();
    Tetrahedral_.clear();
    int i;
    for(i=0;i<nZBins_;i++) {
      sliceSDLists_[i].clear();
    }

    //LOOP OVER ALL FRAMES 
    for (int istep = 0; istep < nFrames; istep += step_) {
      int i;
      for(i=0;i<nZBins_;i++) {
        count_[i]=0;
      } 
        
      reader.readFrame(istep);
      frameCounter_++;
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
      
      if (evaluator_.isDynamic()) {
        seleMan1_.setSelectionSet(evaluator_.evaluate());
        seleMan2_.setSelectionSet(evaluator_.evaluate());
      }
        
      // update the positions of atoms which belong to the rigidbodies
      for (mol = info_->beginMolecule(mi); mol != NULL;
           mol = info_->nextMolecule(mi)) {
        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
             rb = mol->nextRigidBody(rbIter)) {
          rb->updateAtoms();
        }        
      }           

      // outer loop is over the selected StuntDoubles:
      int idk=0;
      for (sd = seleMan1_.beginSelected(isd1); sd != NULL;
           sd = seleMan1_.nextSelected(isd1)) {
        myIndex = sd->getGlobalIndex();
        Qk = 1.0;         
        myNeighbors.clear();
        for(sd2 = seleMan2_.beginSelected(isd2); sd2 != NULL; 
            sd2 = seleMan2_.nextSelected(isd2)){
          if(sd2->getGlobalIndex() != myIndex){
            vec = sd->getPos() - sd2->getPos();       
            if (usePeriodicBoundaryConditions_) 
              currentSnapshot_->wrapVector(vec);
            r = vec.length();             
            
            // Check to see if neighbor is in bond cutoff 
            if (r < rCut_) { 
              myNeighbors.push_back(std::make_pair(r,sd2));
            }
          }
        }
        
        // Sort the vector using predicate and std::sort
        std::sort(myNeighbors.begin(), myNeighbors.end());        
        //std::cerr << myNeighbors.size() <<  " neighbors within " << rCut_  << " A" << " \n";
        // Use only the 4 closest neighbors to do the rest of the work:   
        int nbors =  myNeighbors.size() > 4 ? 4 : myNeighbors.size();
        int nang = int (0.5 * (nbors * (nbors - 1)));
        
        rk = sd->getPos();
        for (int i = 0; i < nbors-1; i++) {
          sdi = myNeighbors[i].second;
          ri = sdi->getPos();
          rik = rk - ri;
          if (usePeriodicBoundaryConditions_) 
            currentSnapshot_->wrapVector(rik);        
          rik.normalize();
          
          for (int j = i+1; j < nbors; j++) {
            sdj = myNeighbors[j].second;
            rj = sdj->getPos();
            rkj = rk - rj;
            if (usePeriodicBoundaryConditions_) 
              currentSnapshot_->wrapVector(rkj);
            rkj.normalize();
            
            cospsi = dot(rik,rkj);
            
            // Calculates scaled Qk for each molecule using calculated
            // angles from 4 or fewer nearest neighbors.
            Qk = Qk - (pow(cospsi + 1.0 / 3.0, 2) * 2.25 / nang);
          }
        }
        
        //std::cerr<<nbors<<endl;
        if (nang > 0) {
          //collectHistogram(Qk);
          
          // Saves positions of StuntDoubles & neighbors with
          // distorted coordination (low Qk value)
          if ((Qk < 0.55) && (Qk > 0.45)) {
            Distorted_.push_back(sd);
            dposition = sd->getPos();
          }
          
          // Saves positions of StuntDoubles & neighbors with
          // tetrahedral coordination (high Qk value)
          if (Qk > 0) {
            Tetrahedral_.push_back(sd);
            tposition = sd->getPos();
          }
          
        }
        
        //wrap the stuntdoubles into a cell      
        Vector3d pos = sd->getPos();
        if (usePeriodicBoundaryConditions_)
          currentSnapshot_->wrapVector(pos);
        sd->setPos(pos);
        // shift molecules by half a box to have bins start at 0
        int binNo = int(nZBins_ * (halfBoxZ_ + pos.z()) / hmat(2,2));
        // Patrick took out the "halfBoxZ_" part in the line above to below
        // int binNo = int(nZBins_ * (pos.z()) / hmat(2,2));
        sliceSDLists_[binNo].push_back(Qk);
        idk++;
      }//outer sd loop
    }//istep loop
 
    //Averaging the value of Qk in each bin
    for(int i=0; i< nZBins_; i++) {
      RealType Qsum=0;
      for (unsigned int k = 0; k < sliceSDLists_[i].size(); ++k) {        
        Qsum=Qsum+sliceSDLists_[i][k];
        count_[i]++;
      }
      //std::cerr<<"past averagin Qk"<<endl;
      //std::cerr<<Qsum<<endl;
      if(count_[i]!=0) {
        Qave_.push_back(Qsum/count_[i]);
      }
      //std::cerr<<count[i]<<endl;
    }
    //std::cerr<<"nZBins_ = "<< nZBins_<<endl;
    //Writing bin#:<Qk> to a file
    std::ofstream rdfStream(outputFilename_.c_str());
    if (rdfStream.is_open()) {
      //rdfStream << "#QkZ\n";
      //rdfStream << "#nFrames:\t" << nProcessed_ << "\n";
      //rdfStream << "#selection: (" << selectionScript_ << ")\n";
      //rdfStream << "#z\tdensity\n";
      for (int i = 0; i < nZBins_; ++i) {
        if(count_[i]!=0) {
          rdfStream << ((hmat(2,2)*i)/nZBins_)+(hmat(2,2)/(2*nZBins_)) 
                    << "\t" << Qave_[i] << "\n";
        }
      }
    }
    
    writeOrderParameter();
    std::cerr << "number of distorted StuntDoubles = " 
              << Distorted_.size() << "\n";
    std::cerr << "number of tetrahedral StuntDoubles = " 
              << Tetrahedral_.size() << "\n";
    collectHistogram(Qk);

  }//void TetrahedralityParam::process() loop
  
  void TetrahedralityParamZ::collectHistogram(RealType Qk) {
    //if (Qk > MinQ_ && Qk < MaxQ_) 
    //  {
    //  int whichBin = int((Qk - MinQ_) / deltaQ_);
    //  Q_histogram_[whichBin] += 1;
    //  }
  }    

  void TetrahedralityParamZ::writeOrderParameter() {
    int nSelected = 0;
    std::ofstream osq((getOutputFileName() + "Q").c_str());
    if (osq.is_open()) {
      osq << "# Tetrahedrality Parameters\n";
      osq << "# selection: (" << selectionScript_ << ")\n";
      osq << "# \n";
      osq.close();
    } else {
      sprintf(painCave.errMsg, "TetrahedralityParamZ: unable to open %s\n", 
              (getOutputFileName() + "q").c_str());
      painCave.isFatal = 1;
      simError();  
    }
    DumpReader reader(info_, dumpFilename_);    
    int nFrames = reader.getNFrames();
    if (nFrames == 1) {
      std::vector<StuntDouble*>::iterator iter;
      std::ofstream osd((getOutputFileName() + "dxyz").c_str());
      if (osd.is_open()) {
        osd << Distorted_.size() << "\n\n";
        
        for (iter = Distorted_.begin(); iter != Distorted_.end(); ++iter) {
          Vector3d position;
          position = (*iter)->getPos();
          osd << "O  " << "\t";
          for (unsigned int z=0; z<position.size(); z++) {
            osd << position[z] << "  " << "\t";
          }
          osd << "\n";
        }
        osd.close();
      }
      std::ofstream ost((getOutputFileName() + "txyz").c_str());
      if (ost.is_open()) {
        ost << Tetrahedral_.size() << "\n\n";       
        for (iter = Tetrahedral_.begin(); iter != Tetrahedral_.end(); ++iter) {
          Vector3d position;            
          position = (*iter)->getPos();
          ost << "O  " << "\t";
          for (unsigned int z=0; z<position.size(); z++) {
            ost << position[z] << "  " << "\t";
          } 
          ost << "\n";
        }
        ost.close();
      }
    }
  }
}


Revision:	1796
Committed:	Mon Sep 10 18:38:44 2012 UTC (12 years, 7 months ago) by gezelter
File size:	12312 byte(s)
Log Message:	Updating MPI calls to C++, fixing a DumpWriter bug, cleaning cruft
#	User	Rev	Content
1	plouden	1762	/*
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, 24107 (2008).
39			* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
41			* Created by J. Daniel Gezelter on 09/26/06.
42			* @author J. Daniel Gezelter
43			* @version $Id: BondOrderParameter.cpp 1442 2010-05-10 17:28:26Z gezelter $
44			*
45			*/
46
47			#include "applications/staticProps/TetrahedralityParamZ.hpp"
48			#include "utils/simError.h"
49			#include "io/DumpReader.hpp"
50			#include "primitives/Molecule.hpp"
51			#include "utils/NumericConstant.hpp"
52			#include <vector>
53			#include <algorithm>
54			#include <fstream>
55
56			using namespace std;
57
58			namespace OpenMD
59			{
60			TetrahedralityParamZ::TetrahedralityParamZ(SimInfo* info,
61			const std::string& filename,
62			const std::string& sele,
63			double rCut, int nzbins) : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan1_(info),seleMan2_(info), nZBins_(nzbins)
64			{
65			//nZBins_ = 50;
66			//std ::cerrnZBins_:"<<nZBins_<<"\t"<<"nzbins:"<<nzbins<<endl;
67			// nZBins_ = 90;
68			//fixed numbe of bins
69	gezelter	1767	count_.resize(nZBins_);
70	plouden	1762	sliceSDLists_.resize(nZBins_);
71			Qave_.resize(nZBins_);
72
73			setOutputName(getPrefix(filename) + ".q");
74
75			evaluator_.loadScriptString(sele);
76			if (!evaluator_.isDynamic())
77			{
78			seleMan1_.setSelectionSet(evaluator_.evaluate());
79			seleMan2_.setSelectionSet(evaluator_.evaluate());
80			}
81
82			// Set up cutoff radius:
83			rCut_ = rCut;
84
85			// Q can take values from 0 to 1
86			MinQ_ = 0.0;
87			MaxQ_ = 1.1;
88			deltaQ_ = (MaxQ_ - MinQ_)/nzbins;
89			}
90
91			TetrahedralityParamZ::~TetrahedralityParamZ()
92			{
93			Q_histogram_.clear();
94			}
95
96			void TetrahedralityParamZ::initalizeHistogram()
97			{
98			std::fill(Q_histogram_.begin(), Q_histogram_.end(), 0);
99			}
100
101
102
103
104			void TetrahedralityParamZ::process()
105			{
106			Molecule* mol;
107			StuntDouble* sd;
108			StuntDouble* sd2;
109			StuntDouble* sdi;
110			StuntDouble* sdj;
111			RigidBody* rb;
112			int myIndex;
113			SimInfo::MoleculeIterator mi;
114			Molecule::RigidBodyIterator rbIter;
115			Vector3d vec;
116			Vector3d ri, rj, rk, rik, rkj, dposition, tposition;
117			RealType r;
118			RealType cospsi;
119			RealType Qk;
120	gezelter	1767
121	plouden	1762	std::vector<std::pair<RealType,StuntDouble*> > myNeighbors;
122			int isd1, isd2;
123			cerr << "After Creation of variables in TP:process()\n";
124			DumpReader reader(info_, dumpFilename_);
125			cerr << "The DumpReader was created?\n";
126			cerr << "nZbins: " << nZBins_ << "\n";
127			int nFrames = reader.getNFrames();
128			frameCounter_ = 0;
129			nProcessed_=nFrames/step_;
130			reader.readFrame(0);
131			currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
132			Mat3x3d hmat = currentSnapshot_->getHmat();
133			zBox_.push_back(hmat(2,2));
134
135			RealType halfBoxZ_ = hmat(2,2) / 2.0;
136
137			Distorted_.clear();
138			Tetrahedral_.clear();
139			int i;
140	gezelter	1796	for(i=0;i<nZBins_;i++) {
141			sliceSDLists_[i].clear();
142			}
143	plouden	1762
144			//LOOP OVER ALL FRAMES
145	gezelter	1796	for (int istep = 0; istep < nFrames; istep += step_) {
146			int i;
147			for(i=0;i<nZBins_;i++) {
148			count_[i]=0;
149			}
150	plouden	1762
151	gezelter	1796	reader.readFrame(istep);
152			frameCounter_++;
153			currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
154
155			if (evaluator_.isDynamic()) {
156			seleMan1_.setSelectionSet(evaluator_.evaluate());
157			seleMan2_.setSelectionSet(evaluator_.evaluate());
158			}
159	plouden	1762
160			// update the positions of atoms which belong to the rigidbodies
161	gezelter	1796	for (mol = info_->beginMolecule(mi); mol != NULL;
162			mol = info_->nextMolecule(mi)) {
163			for (rb = mol->beginRigidBody(rbIter); rb != NULL;
164			rb = mol->nextRigidBody(rbIter)) {
165			rb->updateAtoms();
166			}
167			}
168	plouden	1762
169	gezelter	1796	// outer loop is over the selected StuntDoubles:
170	plouden	1762	int idk=0;
171	gezelter	1796	for (sd = seleMan1_.beginSelected(isd1); sd != NULL;
172			sd = seleMan1_.nextSelected(isd1)) {
173			myIndex = sd->getGlobalIndex();
174			Qk = 1.0;
175			myNeighbors.clear();
176			for(sd2 = seleMan2_.beginSelected(isd2); sd2 != NULL;
177			sd2 = seleMan2_.nextSelected(isd2)){
178			if(sd2->getGlobalIndex() != myIndex){
179			vec = sd->getPos() - sd2->getPos();
180			if (usePeriodicBoundaryConditions_)
181			currentSnapshot_->wrapVector(vec);
182			r = vec.length();
183
184			// Check to see if neighbor is in bond cutoff
185			if (r < rCut_) {
186			myNeighbors.push_back(std::make_pair(r,sd2));
187			}
188			}
189			}
190
191			// Sort the vector using predicate and std::sort
192			std::sort(myNeighbors.begin(), myNeighbors.end());
193			//std::cerr << myNeighbors.size() << " neighbors within " << rCut_ << " A" << " \n";
194			// Use only the 4 closest neighbors to do the rest of the work:
195			int nbors = myNeighbors.size() > 4 ? 4 : myNeighbors.size();
196			int nang = int (0.5 * (nbors * (nbors - 1)));
197
198			rk = sd->getPos();
199			for (int i = 0; i < nbors-1; i++) {
200			sdi = myNeighbors[i].second;
201			ri = sdi->getPos();
202			rik = rk - ri;
203			if (usePeriodicBoundaryConditions_)
204			currentSnapshot_->wrapVector(rik);
205			rik.normalize();
206
207			for (int j = i+1; j < nbors; j++) {
208			sdj = myNeighbors[j].second;
209			rj = sdj->getPos();
210			rkj = rk - rj;
211			if (usePeriodicBoundaryConditions_)
212			currentSnapshot_->wrapVector(rkj);
213			rkj.normalize();
214
215			cospsi = dot(rik,rkj);
216
217			// Calculates scaled Qk for each molecule using calculated
218			// angles from 4 or fewer nearest neighbors.
219			Qk = Qk - (pow(cospsi + 1.0 / 3.0, 2) * 2.25 / nang);
220			}
221			}
222
223			//std::cerr<<nbors<<endl;
224			if (nang > 0) {
225			//collectHistogram(Qk);
226
227			// Saves positions of StuntDoubles & neighbors with
228			// distorted coordination (low Qk value)
229			if ((Qk < 0.55) && (Qk > 0.45)) {
230			Distorted_.push_back(sd);
231			dposition = sd->getPos();
232			}
233
234			// Saves positions of StuntDoubles & neighbors with
235			// tetrahedral coordination (high Qk value)
236			if (Qk > 0) {
237			Tetrahedral_.push_back(sd);
238			tposition = sd->getPos();
239			}
240
241			}
242
243			//wrap the stuntdoubles into a cell
244			Vector3d pos = sd->getPos();
245			if (usePeriodicBoundaryConditions_)
246			currentSnapshot_->wrapVector(pos);
247			sd->setPos(pos);
248			// shift molecules by half a box to have bins start at 0
249			int binNo = int(nZBins_ * (halfBoxZ_ + pos.z()) / hmat(2,2));
250			// Patrick took out the "halfBoxZ_" part in the line above to below
251			// int binNo = int(nZBins_ * (pos.z()) / hmat(2,2));
252			sliceSDLists_[binNo].push_back(Qk);
253			idk++;
254			}//outer sd loop
255			}//istep loop
256	plouden	1762
257			//Averaging the value of Qk in each bin
258	gezelter	1796	for(int i=0; i< nZBins_; i++) {
259			RealType Qsum=0;
260			for (unsigned int k = 0; k < sliceSDLists_[i].size(); ++k) {
261			Qsum=Qsum+sliceSDLists_[i][k];
262			count_[i]++;
263	plouden	1762	}
264	gezelter	1796	//std::cerr<<"past averagin Qk"<<endl;
265			//std::cerr<<Qsum<<endl;
266			if(count_[i]!=0) {
267			Qave_.push_back(Qsum/count_[i]);
268			}
269			//std::cerr<<count[i]<<endl;
270			}
271	plouden	1762	//std::cerr<<"nZBins_ = "<< nZBins_<<endl;
272			//Writing bin#:<Qk> to a file
273			std::ofstream rdfStream(outputFilename_.c_str());
274	gezelter	1796	if (rdfStream.is_open()) {
275			//rdfStream << "#QkZ\n";
276			//rdfStream << "#nFrames:\t" << nProcessed_ << "\n";
277			//rdfStream << "#selection: (" << selectionScript_ << ")\n";
278			//rdfStream << "#z\tdensity\n";
279			for (int i = 0; i < nZBins_; ++i) {
280			if(count_[i]!=0) {
281			rdfStream << ((hmat(2,2)i)/nZBins_)+(hmat(2,2)/(2nZBins_))
282			<< "\t" << Qave_[i] << "\n";
283			}
284	plouden	1762	}
285	gezelter	1796	}
286	plouden	1762
287			writeOrderParameter();
288	gezelter	1796	std::cerr << "number of distorted StuntDoubles = "
289			<< Distorted_.size() << "\n";
290			std::cerr << "number of tetrahedral StuntDoubles = "
291			<< Tetrahedral_.size() << "\n";
292	plouden	1762	collectHistogram(Qk);
293
294			}//void TetrahedralityParam::process() loop
295
296	gezelter	1796	void TetrahedralityParamZ::collectHistogram(RealType Qk) {
297			//if (Qk > MinQ_ && Qk < MaxQ_)
298			// {
299			// int whichBin = int((Qk - MinQ_) / deltaQ_);
300			// Q_histogram_[whichBin] += 1;
301			// }
302	plouden	1762	}
303
304	gezelter	1796	void TetrahedralityParamZ::writeOrderParameter() {
305			int nSelected = 0;
306			std::ofstream osq((getOutputFileName() + "Q").c_str());
307			if (osq.is_open()) {
308			osq << "# Tetrahedrality Parameters\n";
309			osq << "# selection: (" << selectionScript_ << ")\n";
310			osq << "# \n";
311			osq.close();
312			} else {
313			sprintf(painCave.errMsg, "TetrahedralityParamZ: unable to open %s\n",
314			(getOutputFileName() + "q").c_str());
315			painCave.isFatal = 1;
316			simError();
317	plouden	1762	}
318	gezelter	1796	DumpReader reader(info_, dumpFilename_);
319			int nFrames = reader.getNFrames();
320			if (nFrames == 1) {
321			std::vector<StuntDouble*>::iterator iter;
322			std::ofstream osd((getOutputFileName() + "dxyz").c_str());
323			if (osd.is_open()) {
324			osd << Distorted_.size() << "\n\n";
325
326			for (iter = Distorted_.begin(); iter != Distorted_.end(); ++iter) {
327			Vector3d position;
328			position = (*iter)->getPos();
329			osd << "O " << "\t";
330			for (unsigned int z=0; z<position.size(); z++) {
331			osd << position[z] << " " << "\t";
332			}
333			osd << "\n";
334			}
335			osd.close();
336			}
337			std::ofstream ost((getOutputFileName() + "txyz").c_str());
338			if (ost.is_open()) {
339			ost << Tetrahedral_.size() << "\n\n";
340			for (iter = Tetrahedral_.begin(); iter != Tetrahedral_.end(); ++iter) {
341			Vector3d position;
342			position = (*iter)->getPos();
343			ost << "O " << "\t";
344			for (unsigned int z=0; z<position.size(); z++) {
345			ost << position[z] << " " << "\t";
346			}
347			ost << "\n";
348			}
349			ost.close();
350			}
351	plouden	1762	}
352			}
353			}
354
355
356