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

Comparing trunk/src/applications/staticProps/BOPofR.cpp (file contents):
Revision 1785 by jmichalk, Wed Aug 22 18:43:27 2012 UTC vs.
Revision 1992 by gezelter, Thu Apr 24 17:30:00 2014 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 55 | Line 55 | namespace OpenMD {
55   using namespace MATPACK;
56   namespace OpenMD {
57  
58 <  BOPofR::BOPofR(SimInfo* info, const std::string& filename, const std::string& sele, double rCut,
59 <                 int nbins, RealType len) : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info){
58 >  BOPofR::BOPofR(SimInfo* info, const std::string& filename,
59 >                 const std::string& sele, double rCut, int nbins,
60 >                 RealType len) : StaticAnalyser(info, filename),
61 >                                 selectionScript_(sele),
62 >                                 evaluator_(info), seleMan_(info) {
63      
64      setOutputName(getPrefix(filename) + ".bo");
65 <
65 >    
66      evaluator_.loadScriptString(sele);
67      if (!evaluator_.isDynamic()) {
68        seleMan_.setSelectionSet(evaluator_.evaluate());
69      }
70 <
70 >    
71      // Set up cutoff radius and order of the Legendre Polynomial:
72 <
72 >    
73      rCut_ = rCut;
74      nBins_ = nbins;
75      len_ = len;
76 <        
76 >    
77      deltaR_ = len_/nBins_;
78      RCount_.resize(nBins_);
79      WofR_.resize(nBins_);
80      QofR_.resize(nBins_);
81  
82 <        for (int i = 0; i < nBins_; i++){
83 <                RCount_[i] = 0;
84 <                WofR_[i] = 0;
85 <                QofR_[i] = 0;
86 <        }
87 <        
82 >    for (int i = 0; i < nBins_; i++){
83 >      RCount_[i] = 0;
84 >      WofR_[i] = 0;
85 >      QofR_[i] = 0;
86 >    }
87 >    
88      // Make arrays for Wigner3jm
89      RealType* THRCOF = new RealType[2*lMax_+1];
90      // Variables for Wigner routine
# Line 121 | Line 124 | namespace OpenMD {
124    }
125    
126    BOPofR::~BOPofR() {
127 <          /*
128 <        std::cerr << "Freeing stuff" << std::endl;
127 >  /*
128 >    std::cerr << "Freeing stuff" << std::endl;
129      for (int l = 0; l <= lMax_; l++) {
130        for (int m = -l; m <= l; m++) {
131          w3j[std::make_pair(l,m)].clear();
# Line 146 | Line 149 | namespace OpenMD {
149  
150    
151    void BOPofR::initializeHistogram() {
152 <        for (int i = 0; i < nBins_; i++){
153 <                RCount_[i] = 0;
154 <                WofR_[i] = 0;
155 <                QofR_[i] = 0;
156 <        }
152 >    for (int i = 0; i < nBins_; i++){
153 >      RCount_[i] = 0;
154 >      WofR_[i] = 0;
155 >      QofR_[i] = 0;
156 >    }
157    }
158 <
159 <
158 >  
159 >  
160    void BOPofR::process() {
161      Molecule* mol;
162      Atom* atom;
# Line 176 | Line 179 | namespace OpenMD {
179      std::vector<RealType> q2;
180      std::vector<ComplexType> w;
181      std::vector<ComplexType> w_hat;
179 –    std::map<std::pair<int,int>,ComplexType> QBar;
182      std::vector<RealType> Q2;
183      std::vector<RealType> Q;
184      std::vector<ComplexType> W;
185      std::vector<ComplexType> W_hat;
186 <    int nBonds, Nbonds;
186 >    int nBonds;
187      SphericalHarmonic sphericalHarmonic;
188      int i;
189      
# Line 200 | Line 202 | namespace OpenMD {
202      Q.resize(lMax_+1);
203      W.resize(lMax_+1);
204      W_hat.resize(lMax_+1);
203 –    Nbonds = 0;
205  
206      for (int istep = 0; istep < nFrames; istep += step_) {
207        reader.readFrame(istep);
# Line 309 | Line 310 | namespace OpenMD {
310  
311          collectHistogram(q_l, w_hat, distCOM);
312                  
313 < //              printf( "%s  %18.10g %18.10g %18.10g %18.10g \n", sd->getType().c_str(),pos[0],pos[1],pos[2],real(w_hat[6]));
314 <
313 >        //  printf( "%s  %18.10g %18.10g %18.10g %18.10g \n", sd->getType().c_str(),pos[0],pos[1],pos[2],real(w_hat[6]));
314 >        
315        }
316      }
317 <        
317 >    
318      writeOrderParameter();    
319    }
320 +  
321  
320 –  void BOPofR::collectHistogram(std::vector<RealType> q,
321 –                                std::vector<ComplexType> what, RealType distCOM) {
322  
323 +  IcosahedralOfR::IcosahedralOfR(SimInfo* info, const std::string& filename,
324 +                                 const std::string& sele, double rCut,
325 +                                 int nbins, RealType len) : BOPofR(info,
326 +                                                                   filename,
327 +                                                                   sele, rCut,
328 +                                                                   nbins, len) {
329 +  }
330 +
331 +  void IcosahedralOfR::collectHistogram(std::vector<RealType> q,
332 +                                        std::vector<ComplexType> what,
333 +                                        RealType distCOM) {
334 +    
335      if ( distCOM < len_){
336        // Figure out where this distance goes...
337 <      int whichBin = distCOM / deltaR_;
337 >      int whichBin = int(distCOM / deltaR_);
338        RCount_[whichBin]++;
339  
340        if(real(what[6]) < -0.15){                                
341 <                WofR_[whichBin]++;
341 >        WofR_[whichBin]++;
342        }
343 <          if(q[6] > 0.5){
343 >      if(q[6] > 0.5){
344          QofR_[whichBin]++;
345 <          }
346 <    }  
335 <
345 >      }
346 >    }      
347    }
348  
349 <  void BOPofR::writeOrderParameter() {
349 >  FCCOfR::FCCOfR(SimInfo* info, const std::string& filename,
350 >                 const std::string& sele, double rCut,
351 >                 int nbins, RealType len) : BOPofR(info, filename, sele, rCut,
352 >                                                   nbins, len) {}
353 >  
354 >  
355 >  void FCCOfR::collectHistogram(std::vector<RealType> q,
356 >                                std::vector<ComplexType> what,
357 >                                RealType distCOM) {
358      
359 +    if ( distCOM < len_){
360 +      // Figure out where this distance goes...
361 +      int whichBin = int(distCOM / deltaR_);
362 +      RCount_[whichBin]++;
363 +
364 +      if(real(what[4]) < -0.12){                                
365 +        WofR_[whichBin]++;
366 +      }
367 +    }      
368 +  }
369 +  
370 +  void IcosahedralOfR::writeOrderParameter() {
371 +    
372      std::ofstream osq((getOutputFileName() + "qr").c_str());
373  
374      if (osq.is_open()) {
375        
376        // Normalize by number of frames and write it out:
377 <          
377 >      
378        for (int i = 0; i < nBins_; ++i) {
379          RealType Rval = (i + 0.5) * deltaR_;              
380          osq << Rval;
381 <                if (RCount_[i] == 0){
382 <                        osq << "\t" << 0;
383 <                        osq << "\n";
384 <                }else{
385 <                osq << "\t" << (RealType)QofR_[i]/(RealType)RCount_[i];        
386 <                osq << "\n";
387 <                }
381 >        if (RCount_[i] == 0){
382 >          osq << "\t" << 0;
383 >          osq << "\n";
384 >        }else{
385 >          osq << "\t" << (RealType)QofR_[i]/(RealType)RCount_[i];        
386 >          osq << "\n";
387 >        }
388        }
389 <
389 >      
390        osq.close();
391 <
391 >      
392      } else {
393        sprintf(painCave.errMsg, "BOPofR: unable to open %s\n",
394                (getOutputFileName() + "q").c_str());
395        painCave.isFatal = 1;
396        simError();  
397      }
398 <
398 >    
399      std::ofstream osw((getOutputFileName() + "wr").c_str());
400 <
400 >    
401      if (osw.is_open()) {
402        // Normalize by number of frames and write it out:
403        for (int i = 0; i < nBins_; ++i) {
404          RealType Rval = deltaR_ * (i + 0.5);              
405          osw << Rval;
406 <                if (RCount_[i] == 0){
407 <                        osw << "\t" << 0;
408 <                        osw << "\n";
409 <                }else{
410 <                osw << "\t" << (RealType)WofR_[i]/(RealType)RCount_[i];
411 <                osw << "\n";
412 <                }
406 >        if (RCount_[i] == 0){
407 >          osw << "\t" << 0;
408 >          osw << "\n";
409 >        }else{
410 >          osw << "\t" << (RealType)WofR_[i]/(RealType)RCount_[i];
411 >          osw << "\n";
412 >        }
413        }
414 <
414 >      
415        osw.close();
416      } else {
417        sprintf(painCave.errMsg, "BOPofR: unable to open %s\n",
# Line 390 | Line 422 | namespace OpenMD {
422      }
423          
424    }
425 +  void FCCOfR::writeOrderParameter() {
426 +    
427 +    std::ofstream osw((getOutputFileName() + "wr").c_str());
428 +    
429 +    if (osw.is_open()) {
430 +      // Normalize by number of frames and write it out:
431 +      for (int i = 0; i < nBins_; ++i) {
432 +        RealType Rval = deltaR_ * (i + 0.5);              
433 +        osw << Rval;
434 +        if (RCount_[i] == 0){
435 +          osw << "\t" << 0;
436 +          osw << "\n";
437 +        }else{
438 +          osw << "\t" << (RealType)WofR_[i]/(RealType)RCount_[i];
439 +          osw << "\n";
440 +        }
441 +      }
442 +      
443 +      osw.close();
444 +    } else {
445 +      sprintf(painCave.errMsg, "BOPofR: unable to open %s\n",
446 +              (getOutputFileName() + "w").c_str());
447 +      painCave.isFatal = 1;
448 +      simError();  
449 +    
450 +    }
451 +        
452 +  }
453   }

Diff Legend

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