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
    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;
    StuntDouble* sdk;
    RigidBody* rb;
    int myIndex;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
    Molecule::IntegrableObjectIterator ioi;
    Vector3d vec;
    Vector3d ri, rj, rk, rik, rkj, dposition, tposition;
    RealType r;
    RealType dist;
    RealType cospsi;
    RealType Qk;
    RealType Qsum;
    RealType count[nZBins_];
    std::vector<RealType> Qave_;
    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();
          // inner loop is over all StuntDoubles in the system:
          //for (mol = info_->beginMolecule(mi); mol != NULL;mol = info_->nextMolecule(mi)) 
            //{
              //for (sd2 = mol->beginIntegrableObject(ioi); sd2 != NULL; sd2 = mol->nextIntegrableObject(ioi)) 
                //{
                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 (int k = 0; k < sliceSDLists_[i].size(); ++k) 
          {
            Qsum=Qsum+sliceSDLists_[i][k];
            count[i]=count[i]+1;
          }
        //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";
            osd << "1000000.00000000;    34.52893134     0.00000000     0.00000000;     0.00000000    34.52893134     0.00000000;     0.00000000     0.00000000    34.52893134" << "\n";
            
            for (iter = Distorted_.begin(); iter != Distorted_.end(); ++iter) 
              {
                Vector3d position;
                position = (*iter)->getPos();
                osd << "O  " << "\t";
                for (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";
            ost << "1000000.00000000;    34.52893134     0.00000000     0.00000000;     0.00000000    34.52893134     0.00000000;     0.00000000     0.00000000    34.52893134" << "\n";
            
            for (iter = Tetrahedral_.begin(); iter != Tetrahedral_.end(); ++iter) 
              {
                Vector3d position;              
                position = (*iter)->getPos();
                ost << "O  " << "\t";
                for (int z=0; z<position.size(); z++) 
                  {
                    ost << position[z] << "  " << "\t";
                  }
                ost << "\n";
              }
            ost.close();
          }
    }
  }
}


Revision:	1762
Committed:	Thu Jun 28 20:17:33 2012 UTC (12 years, 10 months ago) by plouden
Original Path:	*branches/development/src/applications/staticProps/TetrahedralityParamZ.cpp*
File size:	12388 byte(s)
Log Message:	Calculates the tetrahedrality parameter as a function of zbins. Abandon all hope ye who enter here.
#	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			sliceSDLists_.resize(nZBins_);
70			Qave_.resize(nZBins_);
71
72			setOutputName(getPrefix(filename) + ".q");
73
74			evaluator_.loadScriptString(sele);
75			if (!evaluator_.isDynamic())
76			{
77			seleMan1_.setSelectionSet(evaluator_.evaluate());
78			seleMan2_.setSelectionSet(evaluator_.evaluate());
79			}
80
81			// Set up cutoff radius:
82			rCut_ = rCut;
83
84			// Q can take values from 0 to 1
85			MinQ_ = 0.0;
86			MaxQ_ = 1.1;
87			deltaQ_ = (MaxQ_ - MinQ_)/nzbins;
88			}
89
90			TetrahedralityParamZ::~TetrahedralityParamZ()
91			{
92			Q_histogram_.clear();
93			}
94
95			void TetrahedralityParamZ::initalizeHistogram()
96			{
97			std::fill(Q_histogram_.begin(), Q_histogram_.end(), 0);
98			}
99
100
101
102
103			void TetrahedralityParamZ::process()
104			{
105			Molecule* mol;
106			StuntDouble* sd;
107			StuntDouble* sd2;
108			StuntDouble* sdi;
109			StuntDouble* sdj;
110			StuntDouble* sdk;
111			RigidBody* rb;
112			int myIndex;
113			SimInfo::MoleculeIterator mi;
114			Molecule::RigidBodyIterator rbIter;
115			Molecule::IntegrableObjectIterator ioi;
116			Vector3d vec;
117			Vector3d ri, rj, rk, rik, rkj, dposition, tposition;
118			RealType r;
119			RealType dist;
120			RealType cospsi;
121			RealType Qk;
122			RealType Qsum;
123			RealType count[nZBins_];
124			std::vector<RealType> Qave_;
125			std::vector<std::pair<RealType,StuntDouble*> > myNeighbors;
126			int isd1, isd2;
127			cerr << "After Creation of variables in TP:process()\n";
128			DumpReader reader(info_, dumpFilename_);
129			cerr << "The DumpReader was created?\n";
130			cerr << "nZbins: " << nZBins_ << "\n";
131			int nFrames = reader.getNFrames();
132			frameCounter_ = 0;
133			nProcessed_=nFrames/step_;
134			reader.readFrame(0);
135			currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
136			Mat3x3d hmat = currentSnapshot_->getHmat();
137			zBox_.push_back(hmat(2,2));
138
139			RealType halfBoxZ_ = hmat(2,2) / 2.0;
140
141			Distorted_.clear();
142			Tetrahedral_.clear();
143			int i;
144			for(i=0;i<nZBins_;i++)
145			{
146			sliceSDLists_[i].clear();
147			}
148
149			//LOOP OVER ALL FRAMES
150			for (int istep = 0; istep < nFrames; istep += step_)
151			{
152			int i;
153			for(i=0;i<nZBins_;i++)
154			{
155			count[i]=0;
156			}
157
158			reader.readFrame(istep);
159			frameCounter_++;
160			currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
161
162			if (evaluator_.isDynamic())
163			{
164			seleMan1_.setSelectionSet(evaluator_.evaluate());
165			seleMan2_.setSelectionSet(evaluator_.evaluate());
166			}
167
168			// update the positions of atoms which belong to the rigidbodies
169			for (mol = info_->beginMolecule(mi); mol != NULL;mol = info_->nextMolecule(mi))
170			{
171			for (rb = mol->beginRigidBody(rbIter); rb != NULL;rb = mol->nextRigidBody(rbIter))
172			{
173			rb->updateAtoms();
174			}
175			}
176
177			// outer loop is over the selected StuntDoubles:
178			int idk=0;
179			for (sd = seleMan1_.beginSelected(isd1); sd != NULL;sd = seleMan1_.nextSelected(isd1))
180			{
181			myIndex = sd->getGlobalIndex();
182			Qk = 1.0;
183			myNeighbors.clear();
184			// inner loop is over all StuntDoubles in the system:
185			//for (mol = info_->beginMolecule(mi); mol != NULL;mol = info_->nextMolecule(mi))
186			//{
187			//for (sd2 = mol->beginIntegrableObject(ioi); sd2 != NULL; sd2 = mol->nextIntegrableObject(ioi))
188			//{
189			for(sd2 = seleMan2_.beginSelected(isd2); sd2 != NULL; sd2 = seleMan2_.nextSelected(isd2)){
190			if(sd2->getGlobalIndex() != myIndex){
191			vec = sd->getPos() - sd2->getPos();
192			if (usePeriodicBoundaryConditions_)
193			currentSnapshot_->wrapVector(vec);
194			r = vec.length();
195
196			// Check to see if neighbor is in bond cutoff
197			if (r < rCut_){
198			myNeighbors.push_back(std::make_pair(r,sd2));
199			}
200			}
201			}
202			// }
203			// Sort the vector using predicate and std::sort
204			std::sort(myNeighbors.begin(), myNeighbors.end());
205			//std::cerr << myNeighbors.size() << " neighbors within " << rCut_ << " A" << " \n";
206			// Use only the 4 closest neighbors to do the rest of the work:
207			int nbors = myNeighbors.size()> 4 ? 4 : myNeighbors.size();
208			int nang = int (0.5 * (nbors * (nbors - 1)));
209
210			rk = sd->getPos();
211			for (int i = 0; i < nbors-1; i++)
212			{
213			sdi = myNeighbors[i].second;
214			ri = sdi->getPos();
215			rik = rk - ri;
216			if (usePeriodicBoundaryConditions_)
217			currentSnapshot_->wrapVector(rik);
218			rik.normalize();
219
220			for (int j = i+1; j < nbors; j++)
221			{
222			sdj = myNeighbors[j].second;
223			rj = sdj->getPos();
224			rkj = rk - rj;
225			if (usePeriodicBoundaryConditions_)
226			currentSnapshot_->wrapVector(rkj);
227			rkj.normalize();
228
229			cospsi = dot(rik,rkj);
230
231			// Calculates scaled Qk for each molecule using calculated angles from 4 or fewer nearest neighbors.
232			Qk = Qk - (pow(cospsi + 1.0 / 3.0, 2) * 2.25 / nang);
233			}
234			}
235
236			//std::cerr<<nbors<<endl;
237			if (nang > 0)
238			{
239			//collectHistogram(Qk);
240
241			// Saves positions of StuntDoubles & neighbors with distorted coordination (low Qk value)
242			if ((Qk < 0.55) && (Qk > 0.45))
243			{
244			Distorted_.push_back(sd);
245			dposition = sd->getPos();
246			}
247
248			// Saves positions of StuntDoubles & neighbors with tetrahedral coordination (high Qk value)
249			if (Qk > 0)
250			{
251			Tetrahedral_.push_back(sd);
252			tposition = sd->getPos();
253			}
254
255			}
256
257			//wrap the stuntdoubles into a cell
258			Vector3d pos = sd->getPos();
259			if (usePeriodicBoundaryConditions_)
260			currentSnapshot_->wrapVector(pos);
261			sd->setPos(pos);
262			// shift molecules by half a box to have bins start at 0
263			int binNo = int(nZBins_ * (halfBoxZ_ + pos.z()) / hmat(2,2));
264			//Patrick took out the "halfBoxZ_" part in the line above to below
265			//int binNo = int(nZBins_ * (pos.z()) / hmat(2,2));
266			sliceSDLists_[binNo].push_back(Qk);
267			idk++;
268			}//outer sd loop
269			}//istep loop
270
271			//Averaging the value of Qk in each bin
272			for(int i=0;i< nZBins_; i++)
273			{
274			RealType Qsum=0;
275			for (int k = 0; k < sliceSDLists_[i].size(); ++k)
276			{
277			Qsum=Qsum+sliceSDLists_[i][k];
278			count[i]=count[i]+1;
279			}
280			//std::cerr<<"past averagin Qk"<<endl;
281			//std::cerr<<Qsum<<endl;
282			if(count[i]!=0)
283			{
284			Qave_.push_back(Qsum/count[i]);
285			}
286			//std::cerr<<count[i]<<endl;
287			}
288			//std::cerr<<"nZBins_ = "<< nZBins_<<endl;
289			//Writing bin#:<Qk> to a file
290			std::ofstream rdfStream(outputFilename_.c_str());
291			if (rdfStream.is_open())
292			{
293			//rdfStream << "#QkZ\n";
294			//rdfStream << "#nFrames:\t" << nProcessed_ << "\n";
295			//rdfStream << "#selection: (" << selectionScript_ << ")\n";
296			//rdfStream << "#z\tdensity\n";
297			for (int i = 0; i < nZBins_; ++i)
298			{
299			if(count[i]!=0)
300			{
301			rdfStream << ((hmat(2,2)i)/nZBins_)+(hmat(2,2)/(2nZBins_)) << "\t" << Qave_[i] << "\n";
302			}
303			}
304			}
305
306
307
308
309
310
311			writeOrderParameter();
312			std::cerr << "number of distorted StuntDoubles = " << Distorted_.size() << "\n";
313			std::cerr << "number of tetrahedral StuntDoubles = " << Tetrahedral_.size() << "\n";
314			collectHistogram(Qk);
315
316			}//void TetrahedralityParam::process() loop
317
318			void TetrahedralityParamZ::collectHistogram(RealType Qk)
319			{
320			//if (Qk > MinQ_ && Qk < MaxQ_)
321			// {
322			// int whichBin = int((Qk - MinQ_) / deltaQ_);
323			// Q_histogram_[whichBin] += 1;
324			// }
325			}
326
327			void TetrahedralityParamZ::writeOrderParameter()
328			{
329			int nSelected = 0;
330			std::ofstream osq((getOutputFileName() + "Q").c_str());
331			if (osq.is_open())
332			{
333			osq << "# Tetrahedrality Parameters\n";
334			osq << "# selection: (" << selectionScript_ << ")\n";
335			osq << "# \n";
336			osq.close();
337			}
338			else
339			{
340			sprintf(painCave.errMsg, "TetrahedralityParamZ: unable to open %s\n",
341			(getOutputFileName() + "q").c_str());
342			painCave.isFatal = 1;
343			simError();
344			}
345			DumpReader reader(info_, dumpFilename_);
346			int nFrames = reader.getNFrames();
347			if (nFrames == 1)
348			{
349			std::vector<StuntDouble*>::iterator iter;
350			std::ofstream osd((getOutputFileName() + "dxyz").c_str());
351			if (osd.is_open())
352			{
353			osd << Distorted_.size() << "\n";
354			osd << "1000000.00000000; 34.52893134 0.00000000 0.00000000; 0.00000000 34.52893134 0.00000000; 0.00000000 0.00000000 34.52893134" << "\n";
355
356			for (iter = Distorted_.begin(); iter != Distorted_.end(); ++iter)
357			{
358			Vector3d position;
359			position = (*iter)->getPos();
360			osd << "O " << "\t";
361			for (int z=0; z<position.size(); z++)
362			{
363			osd << position[z] << " " << "\t";
364			}
365			osd << "\n";
366			}
367			osd.close();
368			}
369			std::ofstream ost((getOutputFileName() + "txyz").c_str());
370			if (ost.is_open())
371			{
372			ost << Tetrahedral_.size() << "\n";
373			ost << "1000000.00000000; 34.52893134 0.00000000 0.00000000; 0.00000000 34.52893134 0.00000000; 0.00000000 0.00000000 34.52893134" << "\n";
374
375			for (iter = Tetrahedral_.begin(); iter != Tetrahedral_.end(); ++iter)
376			{
377			Vector3d position;
378			position = (*iter)->getPos();
379			ost << "O " << "\t";
380			for (int z=0; z<position.size(); z++)
381			{
382			ost << position[z] << " " << "\t";
383			}
384			ost << "\n";
385			}
386			ost.close();
387			}
388			}
389			}
390			}
391
392
393