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;
    Molecule::IntegrableObjectIterator ioi;
    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();
          // 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 (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";
            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 (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";
            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 (unsigned int z=0; z<position.size(); z++) 
                  {
                    ost << position[z] << "  " << "\t";
                  }
                ost << "\n";
              }
            ost.close();
          }
    }
  }
}


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