ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/applications/staticProps/BondOrderParameter.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/applications/staticProps/BondOrderParameter.cpp (file contents):
Revision 3007 by gezelter, Tue Sep 19 21:14:11 2006 UTC vs.
Revision 3015 by gezelter, Thu Sep 21 21:47:17 2006 UTC

# Line 51 | Line 51
51   #include "io/DumpReader.hpp"
52   #include "primitives/Molecule.hpp"
53   #include "utils/NumericConstant.hpp"
54 < #include "math/RealSphericalHarmonic.hpp"
54 > #include "math/SphericalHarmonic.hpp"
55 >
56   namespace oopse {
57  
57 –
58    BondOrderParameter::BondOrderParameter(SimInfo* info,
59                                           const std::string& filename,
60                                           const std::string& sele,
61 <                                         double rCut, int lNumber, int nbins)
62 <    : StaticAnalyser(info, filename), selectionScript_(sele),
63 <      evaluator_(info), seleMan_(info){
61 >                                         double rCut, int lNumber, int nbins) : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info){
62      
63 <    setOutputName(getPrefix(filename) + ".obo");
63 >    setOutputName(getPrefix(filename) + ".bo");
64  
65      evaluator_.loadScriptString(sele);
66      if (!evaluator_.isDynamic()) {
# Line 74 | Line 72 | namespace oopse {
72      lNumber_ = lNumber;
73      rCut_ = rCut;
74      mSize_ = 2*lNumber_+1;    
77 –
78 –    // Set the l for the spherical harmonic, it doesn't change
79 –
80 –    sphericalHarmonic.setL(lNumber_);
81 –
82 –    delta_Q = 1.0 / nbins;
83 –    delta_W = 2.0 / nbins;
84 –
85 –    Q_histogram_.resize(nbins);
86 –    W_histogram_.resize(nbins);
87 –
75    }
76  
77 <  void BondOrderParameter::initalizeHistogram() {
91 <    std::fill(Q_histogram_.begin(), Q_histogram_.end(), 0);
92 <    std::fill(W_histogram_.begin(), W_histogram_.end(), 0);
77 >  BondOrderParameter::~BondOrderParameter() {
78    }
79  
80    void BondOrderParameter::process() {
81      Molecule* mol;
82      Atom* atom;
83      RigidBody* rb;
84 +    int myIndex;
85      SimInfo::MoleculeIterator mi;
86      Molecule::RigidBodyIterator rbIter;
87      Molecule::AtomIterator ai;
88      StuntDouble* sd;
89 <    RealType theta;
89 >    Vector3d vec;
90 >    RealType costheta;
91      RealType phi;
92      RealType r;
93      RealType dist;
94 <    std::map<int, RealType> QBar_lm;
94 >    std::map<int,ComplexType> QBar_lm;
95      RealType QSq_l;
96      RealType Q_l;
97 +    ComplexType W_l;
98 +    ComplexType W_l_hat;
99      int nBonds;
100 <    RealSphericalHarmonic sphericalHarmonic;
100 >    SphericalHarmonic sphericalHarmonic;
101      int i, j;
102 <  
103 <  
102 >    // Make arrays for Wigner3jm
103 >    double* THRCOF = new double[mSize_];
104 >    // Variables for Wigner routine
105 >    double l_ = (double)lNumber_;
106 >    double m1Pass, m2Min, m2Max;
107 >    int error, m1, m2, m3;
108 >
109 >    // Set the l for the spherical harmonic, it doesn't change
110 >    sphericalHarmonic.setL(lNumber_);
111 >
112      DumpReader reader(info_, dumpFilename_);    
113      int nFrames = reader.getNFrames();
114 +    frameCounter_ = 0;
115  
118 –
116      for (int istep = 0; istep < nFrames; istep += step_) {
117        reader.readFrame(istep);
118 +      frameCounter_++;
119        currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
120        
121        if (evaluator_.isDynamic()) {
# Line 134 | Line 132 | namespace oopse {
132          }        
133        }      
134        
135 +      nBonds = 0;
136 +      
137 +      for (int m = -lNumber_; m <= lNumber_; m++) {
138 +        QBar_lm[m] = 0.0;
139 +      }
140 +      
141        // outer loop is over the selected StuntDoubles:
142  
143        for (sd = seleMan_.beginSelected(i); sd != NULL;
144             sd = seleMan_.nextSelected(i)) {
145  
146 <        // For this central atom, zero out nBonds and QBar_lm
143 <
144 <        nBonds = 0;
145 <      
146 <        for (int m = -lNumber_; m <= lNumber_; m++) {
147 <          QBar_lm[m] = 0.0;
148 <        }
146 >        myIndex = sd->getGlobalIndex();
147          
148          // inner loop is over all other atoms in the system:
149          
# Line 154 | Line 152 | namespace oopse {
152            for (atom = mol->beginAtom(ai); atom != NULL;
153                 atom = mol->nextAtom(ai)) {
154  
155 +            if (atom->getGlobalIndex() != myIndex) {
156  
157 <            Vector3d vec = sd->getPos() - atom->getPos();      
158 <            currentSnapshot_->wrapVector(vec);
159 <            
160 <            // Calculate "bonds" and build Q_lm(r) where
161 <            //      Q_lm = Y_lm(theta(r),phi(r))                
162 <            // The spherical harmonics are wrt any arbitrary coordinate
163 <            // system, we choose standard spherical coordinates
164 <            
165 <            r = sqrt(pow(vec.x(),2)+pow(vec.y(),2)+pow(vec.z(),2));
166 <            
167 <            // Check to see if neighbor is in bond cutoff
168 <            
169 <            if (r < rCut_) {            
170 <              theta = atan2(vec.y(), vec.x());
171 <              phi = acos(vec.z()/r);
172 <              for(int m = -lNumber_; m <= lNumber_; m++){
173 <                sphericalHarmonic.setM(m);
174 <                QBar_lm[m] += sphericalHarmonic.getValueAt(theta,phi);
175 <              }
176 <              nBonds++;
177 <            }  
178 <          }
179 <        }
181 <        
182 <        // Normalize Qbar
183 <        for (int m = -lNumber_;m <= lNumber_; m++){
184 <          QBar_lm[m] /= nBonds;
185 <        }
186 <
187 <        // Find second order invariant Q_l
188 <
189 <        QSq_l = 0.0;
190 <        for (int m = -lNumber_; m <= lNumber_; m++){
191 <          QSq_l += pow(QBar_lm[m], 2);  
192 <        }
193 <        Q_l = sqrt(QSq_l*(4.0 * NumericConstant::PI / (2.0*(RealType)lNumber_ + 1)));
194 <    
195 <        // Find Third Order Invariant W_l
196 <
197 <        // Make arrays for Wigner3jm
198 <        double* THRCOF = new double[mSize_];
199 <        // Variables for Wigner routine
200 <        double l_ = (double)lNumber_;
201 <        double m2Min, m2Max;
202 <        int error, m1, m2, m3;
203 <        
204 <        W_l_ = 0.0;
205 <        for (int m1 = -lNumber_; m1 <= lNumber_; m1++) {
206 <          // Zero work array
207 <          for (int ii = 0; ii < mSize_; ii+){
208 <            THRCOF[i] = 0.0;
209 <          }
210 <          // Get Wigner coefficients
211 <          Wigner3jm(&l_, &l_, &l_, &(double)m1, &m2Min, &m2Max, THRCOF, &mSize_, &error);
212 <          for (int m_index = 1; i < (int)(m2Max - m2Min-1.0); m_index++) {
213 <            m2 = floor(m2Min) + m_index - 1;
214 <            m3 = -m1-m2;
215 <            W_l_ += THRCOF[m_index]*QBar_lm[m1+lNumber_]*QBar_lm[m2+lNumber_]*QBar_lm[m3+lNumber_];
157 >              vec = sd->getPos() - atom->getPos();      
158 >              currentSnapshot_->wrapVector(vec);
159 >              
160 >              // Calculate "bonds" and build Q_lm(r) where
161 >              //      Q_lm = Y_lm(theta(r),phi(r))                
162 >              // The spherical harmonics are wrt any arbitrary coordinate
163 >              // system, we choose standard spherical coordinates
164 >              
165 >              r = vec.length();
166 >              
167 >              // Check to see if neighbor is in bond cutoff
168 >              
169 >              if (r < rCut_) {
170 >                costheta = vec.z() / r;
171 >                phi = atan2(vec.y(), vec.x());
172 >                
173 >                for(int m = -lNumber_; m <= lNumber_; m++){
174 >                  sphericalHarmonic.setM(m);
175 >                  QBar_lm[m] += sphericalHarmonic.getValueAt(costheta, phi);
176 >                }
177 >                nBonds++;
178 >              }  
179 >            }
180            }
181          }
182 +      }      
183 +    }
184  
185 <        W_l_hat = W_l_ / pow(QSq_l, 1.5);
186 <
187 <        // accumulate histogram data for Q_l and W_l_hat:
222 <
223 <        collectHistogram(Q_l, W_l_hat);
224 <                
225 <      }
185 >    // Normalize Qbar2
186 >    for (int m = -lNumber_;m <= lNumber_; m++){
187 >      QBar_lm[m] /= nBonds;
188      }
189      
190 <    // Normalize by number of frames
190 >    // Find second order invariant Q_l
191 >    
192 >    QSq_l = 0.0;
193      for (int m = -lNumber_; m <= lNumber_; m++){
194 <      QBar_lm[m] /=  nFrames;  
194 >      QSq_l += norm(QBar_lm[m]);
195      }
196      
197 <    
198 <    
235 <    
236 <    
237 <    writeOrderParameter();
197 >    std::cout << "qsl = " << QSq_l << "\n";
198 >    Q_l = sqrt(QSq_l * 4.0 * NumericConstant::PI / (RealType)(2*lNumber_ + 1));
199      
200 <  }
240 <
241 <  
242 <  void BondOrderParameter::processHistogram() {
200 >    // Find Third Order Invariant W_l
201      
202 <    int nPairs = getNPairs();
203 <    RealType volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
204 <    RealType pairDensity = nPairs /volume * 2.0;
205 <    RealType pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0;
206 <
207 <    for(int i = 0 ; i < histogram_.size(); ++i){
208 <
209 <      RealType rLower = i * deltaR_;
210 <      RealType rUpper = rLower + deltaR_;
211 <      RealType volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower );
212 <      RealType nIdeal = volSlice * pairConstant;
213 <
214 <      avgGofr_[i] += histogram_[i] / nIdeal;    
202 >    W_l = 0.0;
203 >    for (int m1 = -lNumber_; m1 <= lNumber_; m1++) {
204 >      // Zero work array
205 >      for (int ii = 0; ii < mSize_; ii++){
206 >        THRCOF[ii] = 0.0;
207 >      }
208 >      // Get Wigner coefficients
209 >      m1Pass = (double)m1;
210 >      
211 >      Wigner3jm(&l_, &l_, &l_,
212 >                &m1Pass, &m2Min, &m2Max,
213 >                THRCOF, &mSize_, &error);
214 >      
215 >      for (int mmm = 0; mmm < (int)(m2Max - m2Min); mmm++) {
216 >        m2 = (int)floor(m2Min) + mmm;
217 >        m3 = -m1-m2;
218 >        W_l += THRCOF[mmm] * QBar_lm[m1] * QBar_lm[m2] * QBar_lm[m3];
219 >      }
220      }
221 <
221 >    
222 >    W_l_hat = W_l / pow(QSq_l, 1.5);              
223 >    
224 >    writeOrderParameter(Q_l, real(W_l_hat));    
225    }
226  
261 –  void BondOrderParameter::collectHistogram(RealType Q_l, RealType W_l_hat) {
227  
228 <    if (Q_l < Max_Q) {
264 <      int whichBin = Q_l / deltaQ_;
265 <      Q_histogram_[whichBin] += 1;
266 <    }
267 <    if (W_l_hat < Max_W) {
268 <      int whichBin = W_l_hat / deltaW_;
269 <      W_histogram_[whichBin] += 1;
270 <    }
271 <  }
272 <  
228 >  void BondOrderParameter::writeOrderParameter(RealType ql, RealType Wlhat) {
229  
274 –  void BondOrderParameter::writeOrderParameter() {
275 –
230      std::ofstream os(getOutputFileName().c_str());
277 –    os << "#Bond Order Parameter\n";
278 –    os << "#selection: (" << selectionScript_ << ")\n";
231  
232 <    for (std::size_t i = 0; i < orderParams_.size(); ++i) {
233 <      os <<  orderParams_[i].p2 << "\t"
234 <         <<  orderParams_[i].director[0] << "\t"
235 <         <<  orderParams_[i].director[1] << "\t"
236 <         <<  orderParams_[i].director[2] << "\t"
237 <         <<  orderParams_[i].angle << "\n";
232 >    if (os.is_open()) {
233 >      
234 >      os << "# Bond Order Parameters\n";
235 >      os << "# selection: (" << selectionScript_ << ")\n";
236 >      os << "# \n";
237 >      os << "# <Q_" << lNumber_ << ">: " << ql << "\n";
238 >      os << "# <W_" << lNumber_ << ">: " << Wlhat << "\n";
239 >      os.close();
240  
241 +    } else {
242 +      sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n",
243 +              getOutputFileName().c_str());
244 +      painCave.isFatal = 1;
245 +      simError();  
246      }
247    }
289 –
290 –
291 –
248   }

Diff Legend

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