ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/staticProps/TetrahedralityParamZ.cpp
(Generate patch)

Comparing trunk/src/applications/staticProps/TetrahedralityParamZ.cpp (file contents):
Revision 1842 by gezelter, Mon Sep 10 20:58:00 2012 UTC vs.
Revision 1843 by gezelter, Tue Jan 29 20:58:08 2013 UTC

# Line 32 | Line 32
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 < *                                                                      
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 < *
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41 > * [6]  Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
42   */
43  
44   #include "applications/staticProps/TetrahedralityParamZ.hpp"
# Line 54 | Line 51 | using namespace std;
51   #include <fstream>
52  
53   using namespace std;
54 <
55 < namespace OpenMD
56 < {
57 <  TetrahedralityParamZ::TetrahedralityParamZ(SimInfo* info,
58 <                                           const std::string& filename,
59 <                                           const std::string& sele,
60 <                                           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 <    count_.resize(nZBins_);
70 <    sliceSDLists_.resize(nZBins_);
71 <    Qave_.resize(nZBins_);
72 <
73 <    setOutputName(getPrefix(filename) + ".q");
54 > namespace OpenMD {
55 >  TetrahedralityParamZ::TetrahedralityParamZ(SimInfo* info,  
56 >                                             const std::string& filename,
57 >                                             const std::string& sele,
58 >                                             double rCut, int nzbins)
59 >    : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info),
60 >      seleMan_(info), nZBins_(nzbins) {
61      
62      evaluator_.loadScriptString(sele);
63 <    if (!evaluator_.isDynamic())
64 <      {
65 <        seleMan1_.setSelectionSet(evaluator_.evaluate());
66 <        seleMan2_.setSelectionSet(evaluator_.evaluate());
67 <      }
81 <
82 <    // Set up cutoff radius:
63 >    if (!evaluator_.isDynamic()) {
64 >      seleMan_.setSelectionSet(evaluator_.evaluate());
65 >    }
66 >    
67 >    // Set up cutoff radius:    
68      rCut_ = rCut;
69  
70 <    // Q can take values from 0 to 1
71 <    MinQ_ = 0.0;
72 <    MaxQ_ = 1.1;
73 <    deltaQ_ = (MaxQ_ - MinQ_)/nzbins;
70 >    // fixed number of bins
71 >    sliceQ_.resize(nZBins_);
72 >    sliceCount_.resize(nZBins_);    
73 >    std::fill(sliceQ_.begin(), sliceQ_.end(), 0.0);
74 >    std::fill(sliceCount_.begin(), sliceCount_.end(), 0);
75 >    
76 >    setOutputName(getPrefix(filename) + ".Qz");
77    }
90 –
91 –  TetrahedralityParamZ::~TetrahedralityParamZ()
92 –  {
93 –    Q_histogram_.clear();
94 –  }
78    
79 <  void TetrahedralityParamZ::initalizeHistogram()
80 <  {
81 <    std::fill(Q_histogram_.begin(), Q_histogram_.end(), 0);
79 >  TetrahedralityParamZ::~TetrahedralityParamZ() {
80 >    sliceQ_.clear();
81 >    sliceCount_.clear();
82 >    zBox_.clear();
83    }
84 <  
85 <
102 <  
103 <
104 <  void TetrahedralityParamZ::process()
105 <  {
84 >    
85 >  void TetrahedralityParamZ::process() {
86      Molecule* mol;
87      StuntDouble* sd;
88      StuntDouble* sd2;
# Line 111 | Line 91 | namespace OpenMD
91      RigidBody* rb;
92      int myIndex;
93      SimInfo::MoleculeIterator mi;
94 +    Molecule::IntegrableObjectIterator ioi;
95      Molecule::RigidBodyIterator rbIter;
96      Vector3d vec;
97 <    Vector3d ri, rj, rk, rik, rkj, dposition, tposition;
97 >    Vector3d ri, rj, rk, rik, rkj;
98      RealType r;
99      RealType cospsi;
100      RealType Qk;
120 –
101      std::vector<std::pair<RealType,StuntDouble*> > myNeighbors;
102 <    int isd1, isd2;
103 <    cerr << "After Creation of variables in TP:process()\n";
102 >    int isd;
103 >
104      DumpReader reader(info_, dumpFilename_);    
125 –    cerr << "The DumpReader was created?\n";
126 –    cerr << "nZbins: " << nZBins_ << "\n";
105      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;
106  
137 –    Distorted_.clear();
138 –    Tetrahedral_.clear();
139 –    int i;
140 –    for(i=0;i<nZBins_;i++) {
141 –      sliceSDLists_[i].clear();
142 –    }
143 –
144 –    //LOOP OVER ALL FRAMES
107      for (int istep = 0; istep < nFrames; istep += step_) {
146 –      int i;
147 –      for(i=0;i<nZBins_;i++) {
148 –        count_[i]=0;
149 –      }
150 –        
108        reader.readFrame(istep);
152 –      frameCounter_++;
109        currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
110        
111 +      Mat3x3d hmat = currentSnapshot_->getHmat();
112 +      zBox_.push_back(hmat(2,2));
113 +      
114 +      RealType halfBoxZ_ = hmat(2,2) / 2.0;      
115 +
116        if (evaluator_.isDynamic()) {
117 <        seleMan1_.setSelectionSet(evaluator_.evaluate());
157 <        seleMan2_.setSelectionSet(evaluator_.evaluate());
117 >        seleMan_.setSelectionSet(evaluator_.evaluate());
118        }
119 <        
119 >      
120        // update the positions of atoms which belong to the rigidbodies
121        for (mol = info_->beginMolecule(mi); mol != NULL;
122             mol = info_->nextMolecule(mi)) {
123          for (rb = mol->beginRigidBody(rbIter); rb != NULL;
124               rb = mol->nextRigidBody(rbIter)) {
125            rb->updateAtoms();
126 <        }        
127 <      }          
128 <
126 >        }
127 >      }
128 >      
129        // outer loop is over the selected StuntDoubles:
130 <      int idk=0;
131 <      for (sd = seleMan1_.beginSelected(isd1); sd != NULL;
132 <           sd = seleMan1_.nextSelected(isd1)) {
130 >      for (sd = seleMan_.beginSelected(isd); sd != NULL;
131 >           sd = seleMan_.nextSelected(isd)) {
132 >        
133          myIndex = sd->getGlobalIndex();
134 +
135          Qk = 1.0;        
136          myNeighbors.clear();
137 <        for(sd2 = seleMan2_.beginSelected(isd2); sd2 != NULL;
138 <            sd2 = seleMan2_.nextSelected(isd2)){
139 <          if(sd2->getGlobalIndex() != myIndex){
140 <            vec = sd->getPos() - sd2->getPos();      
141 <            if (usePeriodicBoundaryConditions_)
142 <              currentSnapshot_->wrapVector(vec);
143 <            r = vec.length();            
137 >        
138 >        // inner loop is over all StuntDoubles in the system:
139 >        
140 >        for (mol = info_->beginMolecule(mi); mol != NULL;
141 >             mol = info_->nextMolecule(mi)) {
142 >
143 >          for (sd2 = mol->beginIntegrableObject(ioi); sd2 != NULL;
144 >               sd2 = mol->nextIntegrableObject(ioi)) {
145              
146 <            // Check to see if neighbor is in bond cutoff
147 <            if (r < rCut_) {
148 <              myNeighbors.push_back(std::make_pair(r,sd2));
146 >            if (sd2->getGlobalIndex() != myIndex) {
147 >              
148 >              vec = sd->getPos() - sd2->getPos();      
149 >              
150 >              if (usePeriodicBoundaryConditions_)
151 >                currentSnapshot_->wrapVector(vec);
152 >              
153 >              r = vec.length();            
154 >
155 >              // Check to see if neighbor is in bond cutoff
156 >              
157 >              if (r < rCut_) {                
158 >                myNeighbors.push_back(std::make_pair(r,sd2));
159 >              }
160              }
161            }
162          }
163          
164          // Sort the vector using predicate and std::sort
165 <        std::sort(myNeighbors.begin(), myNeighbors.end());        
166 <        //std::cerr << myNeighbors.size() <<  " neighbors within " << rCut_  << " A" << " \n";
167 <        // Use only the 4 closest neighbors to do the rest of the work:  
168 <        int nbors =  myNeighbors.size() > 4 ? 4 : myNeighbors.size();
165 >        std::sort(myNeighbors.begin(), myNeighbors.end());
166 >        
167 >        // Use only the 4 closest neighbors to do the rest of the work:
168 >        
169 >        int nbors =  myNeighbors.size()> 4 ? 4 : myNeighbors.size();
170          int nang = int (0.5 * (nbors * (nbors - 1)));
171 <        
171 >        
172          rk = sd->getPos();
173 <        for (int i = 0; i < nbors-1; i++) {
173 >
174 >        for (int i = 0; i < nbors-1; i++) {      
175 >
176            sdi = myNeighbors[i].second;
177            ri = sdi->getPos();
178            rik = rk - ri;
179            if (usePeriodicBoundaryConditions_)
180 <            currentSnapshot_->wrapVector(rik);        
205 <          rik.normalize();
180 >            currentSnapshot_->wrapVector(rik);
181            
182 <          for (int j = i+1; j < nbors; j++) {
182 >          rik.normalize();
183 >
184 >          for (int j = i+1; j < nbors; j++) {      
185 >
186              sdj = myNeighbors[j].second;
187              rj = sdj->getPos();
188              rkj = rk - rj;
# Line 212 | Line 190 | namespace OpenMD
190                currentSnapshot_->wrapVector(rkj);
191              rkj.normalize();
192              
193 <            cospsi = dot(rik,rkj);
194 <            
193 >            cospsi = dot(rik,rkj);          
194 >
195              // Calculates scaled Qk for each molecule using calculated
196              // angles from 4 or fewer nearest neighbors.
197 <            Qk = Qk - (pow(cospsi + 1.0 / 3.0, 2) * 2.25 / nang);
197 >            Qk -=  (pow(cospsi + 1.0 / 3.0, 2) * 2.25 / nang);            
198            }
199          }
200 <        
223 <        //std::cerr<<nbors<<endl;
200 >
201          if (nang > 0) {
202 <          //collectHistogram(Qk);
202 >          if (usePeriodicBoundaryConditions_)
203 >            currentSnapshot_->wrapVector(rk);
204            
205 <          // Saves positions of StuntDoubles & neighbors with
206 <          // distorted coordination (low Qk value)
207 <          if ((Qk < 0.55) && (Qk > 0.45)) {
208 <            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 <
257 <    //Averaging the value of Qk in each bin
258 <    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]++;
205 >          int binNo = int(nZBins_ * (halfBoxZ_ + rk.z()) / hmat(2,2));
206 >          sliceQ_[binNo] += Qk;
207 >          sliceCount_[binNo] += 1;
208 >        }  
209        }
264 –      //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;
210      }
211 <    //std::cerr<<"nZBins_ = "<< nZBins_<<endl;
212 <    //Writing bin#:<Qk> to a file
273 <    std::ofstream rdfStream(outputFilename_.c_str());
274 <    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)/(2*nZBins_))
282 <                    << "\t" << Qave_[i] << "\n";
283 <        }
284 <      }
285 <    }
286 <    
287 <    writeOrderParameter();
288 <    std::cerr << "number of distorted StuntDoubles = "
289 <              << Distorted_.size() << "\n";
290 <    std::cerr << "number of tetrahedral StuntDoubles = "
291 <              << Tetrahedral_.size() << "\n";
292 <    collectHistogram(Qk);
293 <
294 <  }//void TetrahedralityParam::process() loop
211 >    writeQz();
212 >  }
213    
214 <  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 <  }    
214 >  void TetrahedralityParamZ::writeQz() {
215  
216 <  void TetrahedralityParamZ::writeOrderParameter() {
217 <    std::ofstream osq((getOutputFileName() + "Q").c_str());
218 <    if (osq.is_open()) {
219 <      osq << "# Tetrahedrality Parameters\n";
220 <      osq << "# selection: (" << selectionScript_ << ")\n";
221 <      osq << "# \n";
222 <      osq.close();
223 <    } else {
216 >    // compute average box length:
217 >    
218 >    RealType zSum = 0.0;
219 >    for (std::vector<RealType>::iterator j = zBox_.begin();
220 >         j != zBox_.end(); ++j) {
221 >      zSum += *j;      
222 >    }
223 >    RealType zAve = zSum / zBox_.size();
224 >
225 >    std::ofstream qZstream(outputFilename_.c_str());
226 >    if (qZstream.is_open()) {
227 >      qZstream << "#Tetrahedrality Parameters (z)\n";
228 >      qZstream << "#nFrames:\t" << zBox_.size() << "\n";
229 >      qZstream << "#selection: (" << selectionScript_ << ")\n";
230 >      qZstream << "#z\tQk\n";
231 >      for (unsigned int i = 0; i < sliceQ_.size(); ++i) {
232 >        RealType z = zAve * (i+0.5) / sliceQ_.size();
233 >        qZstream << z << "\t" << sliceQ_[i] / sliceCount_[i] << "\n";
234 >      }
235 >      
236 >    } else {      
237        sprintf(painCave.errMsg, "TetrahedralityParamZ: unable to open %s\n",
238 <              (getOutputFileName() + "q").c_str());
238 >              outputFilename_.c_str());
239        painCave.isFatal = 1;
240        simError();  
241 <    }
242 <    DumpReader reader(info_, dumpFilename_);    
318 <    int nFrames = reader.getNFrames();
319 <    if (nFrames == 1) {
320 <      std::vector<StuntDouble*>::iterator iter;
321 <      std::ofstream osd((getOutputFileName() + "dxyz").c_str());
322 <      if (osd.is_open()) {
323 <        osd << Distorted_.size() << "\n\n";
324 <        
325 <        for (iter = Distorted_.begin(); iter != Distorted_.end(); ++iter) {
326 <          Vector3d position;
327 <          position = (*iter)->getPos();
328 <          osd << "O  " << "\t";
329 <          for (unsigned int z=0; z<position.size(); z++) {
330 <            osd << position[z] << "  " << "\t";
331 <          }
332 <          osd << "\n";
333 <        }
334 <        osd.close();
335 <      }
336 <      std::ofstream ost((getOutputFileName() + "txyz").c_str());
337 <      if (ost.is_open()) {
338 <        ost << Tetrahedral_.size() << "\n\n";      
339 <        for (iter = Tetrahedral_.begin(); iter != Tetrahedral_.end(); ++iter) {
340 <          Vector3d position;            
341 <          position = (*iter)->getPos();
342 <          ost << "O  " << "\t";
343 <          for (unsigned int z=0; z<position.size(); z++) {
344 <            ost << position[z] << "  " << "\t";
345 <          }
346 <          ost << "\n";
347 <        }
348 <        ost.close();
349 <      }
350 <    }
241 >    }    
242 >    qZstream.close();
243    }
244   }
245  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines