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 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 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 < *
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41 > * [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)
61 <  {
62 <    //nZBins_ = 50;
63 <    //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& sele1,
58 >                                             const std::string& sele2,
59 >                                             double rCut, int nzbins)
60 >    : StaticAnalyser(info, filename),
61 >      selectionScript1_(sele1), selectionScript2_(sele2),
62 >      seleMan1_(info), seleMan2_(info), evaluator1_(info), evaluator2_(info),
63 >      nZBins_(nzbins) {
64      
65 <    evaluator_.loadScriptString(sele);
66 <    if (!evaluator_.isDynamic())
67 <      {
68 <        seleMan1_.setSelectionSet(evaluator_.evaluate());
69 <        seleMan2_.setSelectionSet(evaluator_.evaluate());
70 <      }
71 <
72 <    // Set up cutoff radius:
65 >    evaluator1_.loadScriptString(sele1);
66 >    if (!evaluator1_.isDynamic()) {
67 >      seleMan1_.setSelectionSet(evaluator1_.evaluate());
68 >    }
69 >    evaluator2_.loadScriptString(sele2);
70 >    if (!evaluator2_.isDynamic()) {
71 >      seleMan2_.setSelectionSet(evaluator2_.evaluate());
72 >    }
73 >    
74 >    // Set up cutoff radius:    
75      rCut_ = rCut;
76  
77 <    // Q can take values from 0 to 1
78 <    MinQ_ = 0.0;
79 <    MaxQ_ = 1.1;
80 <    deltaQ_ = (MaxQ_ - MinQ_)/nzbins;
77 >    // fixed number of bins
78 >    sliceQ_.resize(nZBins_);
79 >    sliceCount_.resize(nZBins_);    
80 >    std::fill(sliceQ_.begin(), sliceQ_.end(), 0.0);
81 >    std::fill(sliceCount_.begin(), sliceCount_.end(), 0);
82 >    
83 >    setOutputName(getPrefix(filename) + ".Qz");
84    }
90 –
91 –  TetrahedralityParamZ::~TetrahedralityParamZ()
92 –  {
93 –    Q_histogram_.clear();
94 –  }
85    
86 <  void TetrahedralityParamZ::initalizeHistogram()
87 <  {
88 <    std::fill(Q_histogram_.begin(), Q_histogram_.end(), 0);
86 >  TetrahedralityParamZ::~TetrahedralityParamZ() {
87 >    sliceQ_.clear();
88 >    sliceCount_.clear();
89 >    zBox_.clear();
90    }
91 <  
92 <
102 <  
103 <
104 <  void TetrahedralityParamZ::process()
105 <  {
91 >    
92 >  void TetrahedralityParamZ::process() {
93      Molecule* mol;
94      StuntDouble* sd;
95      StuntDouble* sd2;
# Line 112 | Line 99 | namespace OpenMD
99      int myIndex;
100      SimInfo::MoleculeIterator mi;
101      Molecule::RigidBodyIterator rbIter;
115 –    Molecule::IntegrableObjectIterator ioi;
102      Vector3d vec;
103 <    Vector3d ri, rj, rk, rik, rkj, dposition, tposition;
103 >    Vector3d ri, rj, rk, rik, rkj;
104      RealType r;
105      RealType cospsi;
106      RealType Qk;
121 –
107      std::vector<std::pair<RealType,StuntDouble*> > myNeighbors;
108 <    int isd1, isd2;
109 <    cerr << "After Creation of variables in TP:process()\n";
108 >    int isd1;
109 >    int isd2;
110 >
111      DumpReader reader(info_, dumpFilename_);    
126 –    cerr << "The DumpReader was created?\n";
127 –    cerr << "nZbins: " << nZBins_ << "\n";
112      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;
113  
114 <    Distorted_.clear();
115 <    Tetrahedral_.clear();
116 <    int i;
117 <    for(i=0;i<nZBins_;i++)
118 <      {
119 <        sliceSDLists_[i].clear();
120 <      }
114 >    for (int istep = 0; istep < nFrames; istep += step_) {
115 >      reader.readFrame(istep);
116 >      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
117 >      
118 >      Mat3x3d hmat = currentSnapshot_->getHmat();
119 >      zBox_.push_back(hmat(2,2));
120 >      
121 >      RealType halfBoxZ_ = hmat(2,2) / 2.0;      
122  
123 <    //LOOP OVER ALL FRAMES
124 <    for (int istep = 0; istep < nFrames; istep += step_)
125 <      {
126 <        int i;
127 <        for(i=0;i<nZBins_;i++)
128 <          {
129 <            count_[i]=0;
130 <          }    
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 <        
123 >      if (evaluator1_.isDynamic()) {
124 >        seleMan1_.setSelectionSet(evaluator1_.evaluate());
125 >      }
126 >      
127 >      if (evaluator2_.isDynamic()) {
128 >        seleMan2_.setSelectionSet(evaluator2_.evaluate());
129 >      }
130 >      
131        // update the positions of atoms which belong to the rigidbodies
132 <      for (mol = info_->beginMolecule(mi); mol != NULL;mol = info_->nextMolecule(mi))
133 <        {
134 <          for (rb = mol->beginRigidBody(rbIter); rb != NULL;rb = mol->nextRigidBody(rbIter))
135 <            {
136 <              rb->updateAtoms();
137 <            }        
138 <        }          
132 >      for (mol = info_->beginMolecule(mi); mol != NULL;
133 >           mol = info_->nextMolecule(mi)) {
134 >        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
135 >             rb = mol->nextRigidBody(rbIter)) {
136 >          rb->updateAtoms();
137 >        }
138 >      }
139 >      
140 >      // outer loop is over the selected StuntDoubles:
141 >      for (sd = seleMan1_.beginSelected(isd1); sd != NULL;
142 >           sd = seleMan1_.nextSelected(isd1)) {
143 >        
144 >        myIndex = sd->getGlobalIndex();
145 >        
146 >        Qk = 1.0;        
147 >        myNeighbors.clear();      
148  
149 <       // outer loop is over the selected StuntDoubles:
150 <      int idk=0;
151 <      for (sd = seleMan1_.beginSelected(isd1); sd != NULL;sd = seleMan1_.nextSelected(isd1))
152 <        {
153 <          myIndex = sd->getGlobalIndex();
154 <          Qk = 1.0;      
155 <          myNeighbors.clear();
156 <          // inner loop is over all StuntDoubles in the system:
157 <          //for (mol = info_->beginMolecule(mi); mol != NULL;mol = info_->nextMolecule(mi))
158 <            //{
159 <              //for (sd2 = mol->beginIntegrableObject(ioi); sd2 != NULL; sd2 = mol->nextIntegrableObject(ioi))
160 <                //{
161 <                for(sd2 = seleMan2_.beginSelected(isd2); sd2 != NULL; sd2 = seleMan2_.nextSelected(isd2)){
162 <                        if(sd2->getGlobalIndex() != myIndex){
163 <                        vec = sd->getPos() - sd2->getPos();      
164 <                        if (usePeriodicBoundaryConditions_)
165 <                                        currentSnapshot_->wrapVector(vec);
166 <                        r = vec.length();            
167 <
168 <                        // Check to see if neighbor is in bond cutoff
169 <                        if (r < rCut_){
170 <                                        myNeighbors.push_back(std::make_pair(r,sd2));
171 <                                }
172 <                        }
173 <                }
174 <           // }
175 <          // Sort the vector using predicate and std::sort
176 <          std::sort(myNeighbors.begin(), myNeighbors.end());      
177 <          //std::cerr << myNeighbors.size() <<  " neighbors within " << rCut_  << " A" << " \n";
178 <          // Use only the 4 closest neighbors to do the rest of the work:        
179 <          int nbors =  myNeighbors.size()> 4 ? 4 : myNeighbors.size();
180 <          int nang = int (0.5 * (nbors * (nbors - 1)));
181 <          
182 <          rk = sd->getPos();
183 <          for (int i = 0; i < nbors-1; i++)
184 <            {                
185 <              sdi = myNeighbors[i].second;
186 <              ri = sdi->getPos();
187 <              rik = rk - ri;
188 <              if (usePeriodicBoundaryConditions_)
189 <                currentSnapshot_->wrapVector(rik);            
190 <              rik.normalize();
191 <              
192 <              for (int j = i+1; j < nbors; j++)
193 <                {                    
194 <                  sdj = myNeighbors[j].second;
195 <                  rj = sdj->getPos();
196 <                  rkj = rk - rj;
197 <                  if (usePeriodicBoundaryConditions_)
198 <                    currentSnapshot_->wrapVector(rkj);
199 <                  rkj.normalize();
200 <                  
201 <                  cospsi = dot(rik,rkj);
202 <                  
203 <                  // Calculates scaled Qk for each molecule using calculated angles from 4 or fewer nearest neighbors.
204 <                  Qk = Qk - (pow(cospsi + 1.0 / 3.0, 2) * 2.25 / nang);
205 <                }
206 <            }
207 <          
208 <          //std::cerr<<nbors<<endl;
209 <           if (nang > 0)
210 <           {
211 <             //collectHistogram(Qk);
212 <              
213 <              // 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 <        for (unsigned int k = 0; k < sliceSDLists_[i].size(); ++k)
273 <          {
274 <            Qsum=Qsum+sliceSDLists_[i][k];
275 <            count_[i]++;
276 <          }
277 <        //std::cerr<<"past averagin Qk"<<endl;
278 <        //std::cerr<<Qsum<<endl;
279 <        if(count_[i]!=0)
280 <          {
281 <            Qave_.push_back(Qsum/count_[i]);
282 <          }
283 <        //std::cerr<<count[i]<<endl;
149 >        for (sd2 = seleMan2_.beginSelected(isd2); sd2 != NULL;
150 >             sd2 = seleMan2_.nextSelected(isd2)) {
151 >          
152 >          if (sd2->getGlobalIndex() != myIndex) {
153 >            
154 >            vec = sd->getPos() - sd2->getPos();      
155 >            
156 >            if (usePeriodicBoundaryConditions_)
157 >              currentSnapshot_->wrapVector(vec);
158 >            
159 >            r = vec.length();            
160 >            
161 >            // Check to see if neighbor is in bond cutoff
162 >            
163 >            if (r < rCut_) {                
164 >              myNeighbors.push_back(std::make_pair(r,sd2));
165 >            }
166 >          }
167 >        }
168 >        
169 >        // Sort the vector using predicate and std::sort
170 >        std::sort(myNeighbors.begin(), myNeighbors.end());
171 >        
172 >        // Use only the 4 closest neighbors to do the rest of the work:
173 >        
174 >        int nbors =  myNeighbors.size()> 4 ? 4 : myNeighbors.size();
175 >        int nang = int (0.5 * (nbors * (nbors - 1)));
176 >        
177 >        rk = sd->getPos();
178 >        
179 >        for (int i = 0; i < nbors-1; i++) {      
180 >          
181 >          sdi = myNeighbors[i].second;
182 >          ri = sdi->getPos();
183 >          rik = rk - ri;
184 >          if (usePeriodicBoundaryConditions_)
185 >            currentSnapshot_->wrapVector(rik);
186 >          
187 >          rik.normalize();
188 >          
189 >          for (int j = i+1; j < nbors; j++) {      
190 >            
191 >            sdj = myNeighbors[j].second;
192 >            rj = sdj->getPos();
193 >            rkj = rk - rj;
194 >            if (usePeriodicBoundaryConditions_)
195 >              currentSnapshot_->wrapVector(rkj);
196 >            rkj.normalize();
197 >            
198 >            cospsi = dot(rik,rkj);          
199 >            
200 >            // Calculates scaled Qk for each molecule using calculated
201 >            // angles from 4 or fewer nearest neighbors.
202 >            Qk -=  (pow(cospsi + 1.0 / 3.0, 2) * 2.25 / nang);            
203 >          }
204 >        }
205 >        
206 >        if (nang > 0) {
207 >          if (usePeriodicBoundaryConditions_)
208 >            currentSnapshot_->wrapVector(rk);
209 >          
210 >          int binNo = int(nZBins_ * (halfBoxZ_ + rk.z()) / hmat(2,2));
211 >          sliceQ_[binNo] += Qk;
212 >          sliceCount_[binNo] += 1;
213 >        }  
214        }
215 <    //std::cerr<<"nZBins_ = "<< nZBins_<<endl;
216 <    //Writing bin#:<Qk> to a file
217 <    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 <            if(count_[i]!=0)
297 <              {
298 <                rdfStream << ((hmat(2,2)*i)/nZBins_)+(hmat(2,2)/(2*nZBins_)) << "\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
215 >    }
216 >    writeQz();
217 >  }
218    
219 <  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 <  }    
219 >  void TetrahedralityParamZ::writeQz() {
220  
221 <  void TetrahedralityParamZ::writeOrderParameter()
222 <  {  
223 <   int nSelected = 0;
224 <  std::ofstream osq((getOutputFileName() + "Q").c_str());
225 <  if (osq.is_open())
226 <    {
330 <        osq << "# Tetrahedrality Parameters\n";
331 <        osq << "# selection: (" << selectionScript_ << ")\n";
332 <        osq << "# \n";
333 <        osq.close();
221 >    // compute average box length:
222 >    
223 >    RealType zSum = 0.0;
224 >    for (std::vector<RealType>::iterator j = zBox_.begin();
225 >         j != zBox_.end(); ++j) {
226 >      zSum += *j;      
227      }
228 <  else
229 <    {
230 <        sprintf(painCave.errMsg, "TetrahedralityParamZ: unable to open %s\n",
231 <                (getOutputFileName() + "q").c_str());
232 <        painCave.isFatal = 1;
233 <        simError();  
234 <    }
235 <  DumpReader reader(info_, dumpFilename_);    
236 <  int nFrames = reader.getNFrames();
237 <  if (nFrames == 1)
238 <    {
239 <        std::vector<StuntDouble*>::iterator iter;
240 <        std::ofstream osd((getOutputFileName() + "dxyz").c_str());
241 <        if (osd.is_open())
242 <          {
243 <            osd << Distorted_.size() << "\n";
244 <            osd << "1000000.00000000;    34.52893134     0.00000000     0.00000000;     0.00000000    34.52893134     0.00000000;     0.00000000     0.00000000    34.52893134" << "\n";
245 <            
246 <            for (iter = Distorted_.begin(); iter != Distorted_.end(); ++iter)
247 <              {
248 <                Vector3d position;
249 <                position = (*iter)->getPos();
250 <                osd << "O  " << "\t";
358 <                for (unsigned int z=0; z<position.size(); z++)
359 <                  {
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 <                for (unsigned int z=0; z<position.size(); z++)
378 <                  {
379 <                    ost << position[z] << "  " << "\t";
380 <                  }
381 <                ost << "\n";
382 <              }
383 <            ost.close();
384 <          }
385 <    }
228 >    RealType zAve = zSum / zBox_.size();
229 >
230 >    std::ofstream qZstream(outputFilename_.c_str());
231 >    if (qZstream.is_open()) {
232 >      qZstream << "#Tetrahedrality Parameters (z)\n";
233 >      qZstream << "#nFrames:\t" << zBox_.size() << "\n";
234 >      qZstream << "#selection 1: (" << selectionScript1_ << ")\n";
235 >      qZstream << "#selection 2: (" << selectionScript2_ << ")\n";
236 >      qZstream << "#z\tQk\n";
237 >      for (unsigned int i = 0; i < sliceQ_.size(); ++i) {
238 >        RealType z = zAve * (i+0.5) / sliceQ_.size();
239 >        if (sliceCount_[i] != 0) {
240 >          qZstream << z << "\t" << sliceQ_[i] / sliceCount_[i] << "\n";
241 >        }
242 >      }
243 >      
244 >    } else {      
245 >      sprintf(painCave.errMsg, "TetrahedralityParamZ: unable to open %s\n",
246 >              outputFilename_.c_str());
247 >      painCave.isFatal = 1;
248 >      simError();  
249 >    }    
250 >    qZstream.close();
251    }
252   }
253  

Diff Legend

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