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

Comparing trunk/src/applications/staticProps/TetrahedralityParam.cpp (file contents):
Revision 1524 by kstocke1, Fri Nov 19 20:48:18 2010 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

# Line 35 | Line 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]  Vardeman & Gezelter, in progress (2009).                        
40 < *
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (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 $
# Line 56 | Line 56 | namespace OpenMD {
56    TetrahedralityParam::TetrahedralityParam(SimInfo* info,
57                                             const std::string& filename,
58                                             const std::string& sele,
59 <                                           double rCut, int nbins) : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info){
59 >                                           double rCut, int nbins) :
60 >    StaticAnalyser(info, filename), selectionScript_(sele),
61 >    seleMan_(info), evaluator_(info) {
62      
63      setOutputName(getPrefix(filename) + ".q");
64 <
64 >    
65      evaluator_.loadScriptString(sele);
66      if (!evaluator_.isDynamic()) {
67        seleMan_.setSelectionSet(evaluator_.evaluate());
# Line 84 | Line 86 | namespace OpenMD {
86      Q_histogram_.clear();
87    }
88    
89 <  void TetrahedralityParam::initalizeHistogram() {
89 >  void TetrahedralityParam::initializeHistogram() {
90      std::fill(Q_histogram_.begin(), Q_histogram_.end(), 0);
91    }
92    
# Line 94 | Line 96 | namespace OpenMD {
96      StuntDouble* sd2;
97      StuntDouble* sdi;
98      StuntDouble* sdj;
97 –    StuntDouble* sdk;
99      RigidBody* rb;
100      int myIndex;
101      SimInfo::MoleculeIterator mi;
# Line 103 | Line 104 | namespace OpenMD {
104      Vector3d vec;
105      Vector3d ri, rj, rk, rik, rkj, dposition, tposition;
106      RealType r;
106 –    RealType dist;
107      RealType cospsi;
108      RealType Qk;
109      std::vector<std::pair<RealType,StuntDouble*> > myNeighbors;
# Line 136 | Line 136 | namespace OpenMD {
136        }          
137              
138  
139 <       // outer loop is over the selected StuntDoubles:
139 >      // outer loop is over the selected StuntDoubles:
140  
141        for (sd = seleMan_.beginSelected(isd); sd != NULL;
142             sd = seleMan_.nextSelected(isd)) {
# Line 176 | Line 176 | namespace OpenMD {
176          // Sort the vector using predicate and std::sort
177          std::sort(myNeighbors.begin(), myNeighbors.end());
178  
179 <        std::cerr << myNeighbors.size() <<  " neighbors within " << rCut_  << " A" << " \n";
179 >        //std::cerr << myNeighbors.size() <<  " neighbors within "
180 >        //          << rCut_  << " A" << " \n";
181          
182          // Use only the 4 closest neighbors to do the rest of the work:
183          
184 <        int nbors =  myNeighbors.size();
184 <        // > 4 ? 4 : myNeighbors.size();
184 >        int nbors =  myNeighbors.size()> 4 ? 4 : myNeighbors.size();
185          int nang = int (0.5 * (nbors * (nbors - 1)));
186  
187          rk = sd->getPos();
188 <
188 >        //std::cerr<<nbors<<endl;
189          for (int i = 0; i < nbors-1; i++) {      
190  
191            sdi = myNeighbors[i].second;
# Line 209 | Line 209 | namespace OpenMD {
209  
210              //std::cerr << "cos(psi) = " << cospsi << " \n";
211  
212 <            // Calculates scaled Qk for each molecule using calculated angles from 4 or fewer nearest neighbors.
212 >            // Calculates scaled Qk for each molecule using calculated
213 >            // angles from 4 or fewer nearest neighbors.
214              Qk = Qk - (pow(cospsi + 1.0 / 3.0, 2) * 2.25 / nang);
215 <          
215 >            //std::cerr<<Qk<<"\t"<<nang<<endl;
216            }
217          }
218 <
218 >        //std::cerr<<nang<<endl;
219          if (nang > 0) {
220            collectHistogram(Qk);
221  
222 <        // Saves positions of StuntDoubles & neighbors with distorted coordination (low Qk value)
223 <        if ((Qk < 0.55) && (Qk > 0.45)) {
224 <
225 <          Distorted_.push_back(sd);
226 <
227 <          dposition = sd->getPos();
228 <          //std::cerr << "distorted position \t" << dposition << "\n";
229 <        }
230 <
230 <        // Saves positions of StuntDoubles & neighbors with tetrahedral coordination (high Qk value)
231 <        if (Qk > 0.95) {
222 >          // Saves positions of StuntDoubles & neighbors with distorted
223 >          // coordination (low Qk value)
224 >          if ((Qk < 0.55) && (Qk > 0.45)) {
225 >            //std::cerr<<Distorted_.size()<<endl;
226 >            Distorted_.push_back(sd);
227 >            //std::cerr<<Distorted_.size()<<endl;
228 >            dposition = sd->getPos();
229 >            //std::cerr << "distorted position \t" << dposition << "\n";
230 >          }
231  
232 <          Tetrahedral_.push_back(sd);
232 >          // Saves positions of StuntDoubles & neighbors with
233 >          // tetrahedral coordination (high Qk value)
234 >          if (Qk > 0.05) {
235  
236 <          tposition = sd->getPos();
236 <          //std::cerr << "tetrahedral position \t" << tposition << "\n";
237 <        }
236 >            Tetrahedral_.push_back(sd);
237  
238 +            tposition = sd->getPos();
239 +            //std::cerr << "tetrahedral position \t" << tposition << "\n";
240 +          }
241 +
242 +          //std::cerr<<Tetrahedral_.size()<<endl;
243 +      
244          }
245  
246        }
247      }
248      
249      writeOrderParameter();
250 <    std::cerr << "number of distorted StuntDoubles = " << Distorted_.size() << "\n";
251 <    std::cerr << "number of tetrahedral StuntDoubles = " << Tetrahedral_.size() << "\n";
250 >    std::cerr << "number of distorted StuntDoubles = "
251 >              << Distorted_.size() << "\n";
252 >    std::cerr << "number of tetrahedral StuntDoubles = "
253 >              << Tetrahedral_.size() << "\n";
254    }
255          
256    void TetrahedralityParam::collectHistogram(RealType Qk) {
# Line 261 | Line 268 | namespace OpenMD {
268      int nSelected = 0;
269  
270      for (int i = 0; i < nBins_; ++i) {
271 <      nSelected = nSelected + Q_histogram_[i]*deltaQ_;
271 >      nSelected = nSelected + int(Q_histogram_[i]*deltaQ_);
272      }
273 <
273 >    
274      std::ofstream osq((getOutputFileName() + "Q").c_str());
275  
276      if (osq.is_open()) {
# Line 293 | Line 300 | namespace OpenMD {
300  
301      if (nFrames == 1) {
302  
303 <    std::vector<StuntDouble*>::iterator iter;
304 <    std::ofstream osd((getOutputFileName() + "dxyz").c_str());
298 <
299 <    if (osd.is_open()) {
303 >      std::vector<StuntDouble*>::iterator iter;
304 >      std::ofstream osd((getOutputFileName() + "dxyz").c_str());
305  
306 <      osd << Distorted_.size() << "\n";
306 >      if (osd.is_open()) {
307  
308 <      osd << "1000000.00000000;    34.52893134     0.00000000     0.00000000;     0.00000000    34.52893134     0.00000000;     0.00000000     0.00000000    34.52893134" << "\n";
308 >        osd << Distorted_.size() << "\n";
309 >        osd << "\n";
310        
311 <      for (iter = Distorted_.begin(); iter != Distorted_.end(); ++iter) {
312 <
313 <        Vector3d position;
314 <
315 <        position = (*iter)->getPos();
310 <
311 <        osd << "O  " << "\t";
312 <
313 <          for (int z=0; z<position.size(); z++) {
314 <
311 >        for (iter = Distorted_.begin(); iter != Distorted_.end(); ++iter) {
312 >          Vector3d position;
313 >          position = (*iter)->getPos();
314 >          osd << "O  " << "\t";
315 >          for (unsigned int z = 0; z < position.size(); z++) {
316              osd << position[z] << "  " << "\t";
317            }
317 –
318            osd << "\n";
319 <
319 >        }
320 >        osd.close();
321        }
322  
322 –      osd.close();
323 –    }
323  
324 <
326 <    std::ofstream ost((getOutputFileName() + "txyz").c_str());
324 >      std::ofstream ost((getOutputFileName() + "txyz").c_str());
325      
326 <    if (ost.is_open()) {
326 >      if (ost.is_open()) {
327  
328 <      ost << Tetrahedral_.size() << "\n";
329 <
332 <      ost << "1000000.00000000;    34.52893134     0.00000000     0.00000000;     0.00000000    34.52893134     0.00000000;     0.00000000     0.00000000    34.52893134" << "\n";
328 >        ost << Tetrahedral_.size() << "\n";
329 >        ost << "\n";
330        
331 <      for (iter = Tetrahedral_.begin(); iter != Tetrahedral_.end(); ++iter) {
331 >        for (iter = Tetrahedral_.begin(); iter != Tetrahedral_.end(); ++iter) {
332 >          Vector3d position;
333 >          position = (*iter)->getPos();
334  
335 <        Vector3d position;
335 >          ost << "O  " << "\t";
336  
337 <        position = (*iter)->getPos();
339 <
340 <        ost << "O  " << "\t";
341 <
342 <          for (int z=0; z<position.size(); z++) {
343 <
337 >          for (unsigned int z = 0; z < position.size(); z++) {
338              ost << position[z] << "  " << "\t";
339            }
346 –
340            ost << "\n";
341 <
341 >        }
342 >        ost.close();
343        }
350 –
351 –      ost.close();
344      }
353 –    
345    }
346   }
356 – }
357 –      
347  
348  
349 +

Diff Legend

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