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 1785 by jmichalk, Wed Aug 22 18:43:27 2012 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).          
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.
# 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 134 | 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 174 | 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          
# Line 182 | Line 185 | namespace OpenMD {
185          int nang = int (0.5 * (nbors * (nbors - 1)));
186  
187          rk = sd->getPos();
188 <        std::cerr<<nbors<<endl;
188 >        //std::cerr<<nbors<<endl;
189          for (int i = 0; i < nbors-1; i++) {      
190  
191            sdi = myNeighbors[i].second;
# Line 206 | 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              //std::cerr<<Qk<<"\t"<<nang<<endl;
216            }
# Line 215 | Line 219 | namespace OpenMD {
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 <          //std::cerr<<Distorted_.size()<<endl;
225 <          Distorted_.push_back(sd);
226 <          //std::cerr<<Distorted_.size()<<endl;
227 <          dposition = sd->getPos();
228 <          //std::cerr << "distorted position \t" << dposition << "\n";
229 <        }
230 <
227 <        // Saves positions of StuntDoubles & neighbors with tetrahedral coordination (high Qk value)
228 <        if (Qk > 0.05) {
229 <
230 <          Tetrahedral_.push_back(sd);
231 <
232 <          tposition = sd->getPos();
233 <          //std::cerr << "tetrahedral position \t" << tposition << "\n";
234 <        }
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 <        //std::cerr<<Tetrahedral_.size()<<endl;
232 >          // Saves positions of StuntDoubles & neighbors with
233 >          // tetrahedral coordination (high Qk value)
234 >          if (Qk > 0.05) {
235  
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      
274      std::ofstream osq((getOutputFileName() + "Q").c_str());
# 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 <        Vector3d position;
313 <        position = (*iter)->getPos();
314 <        osd << "O  " << "\t";
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            }
318            osd << "\n";
319 +        }
320 +        osd.close();
321        }
314 –      osd.close();
315 –    }
322  
323  
324 <    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 <
324 <      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  
330 –        position = (*iter)->getPos();
331 –
332 –        ost << "O  " << "\t";
333 –
337            for (unsigned int z = 0; z < position.size(); z++) {
335 –
338              ost << position[z] << "  " << "\t";
339            }
338 –
340            ost << "\n";
341 <
341 >        }
342 >        ost.close();
343        }
342 –
343 –      ost.close();
344      }
345 –    
345    }
346   }
348 – }
349 –      
347  
348  
349 +

Diff Legend

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