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

Comparing trunk/src/applications/staticProps/BondOrderParameter.cpp (file contents):
Revision 1038 by chuckv, Mon Sep 18 21:31:23 2006 UTC vs.
Revision 1039 by gezelter, Tue Sep 19 21:14:11 2006 UTC

# Line 45 | Line 45
45   *     Phys Rev B, 28,784,1983
46   *
47   */
48 <
48 >
49   #include "applications/staticProps/BondOrderParameter.hpp"
50   #include "utils/simError.h"
51   #include "io/DumpReader.hpp"
# Line 55 | Line 55 | namespace oopse {
55   namespace oopse {
56  
57  
58 <  BondOrderParameter::BondOrderParameter(SimInfo* info, const std::string& filename, const std::string& sele1,
59 <                                         const std::string& sele2, double rCut, int lNumber)
60 <    : StaticAnalyser(info, filename),
61 <      selectionScript1_(sele1), evaluator1_(info),
62 <      seleMan1_(info){
63 <
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){
64 >    
65      setOutputName(getPrefix(filename) + ".obo");
66  
67 <    evaluator1_.loadScriptString(sele1);
68 <    evaluator2_.loadScriptString(sele2);
69 <
69 <    if (!evaluator1_.isDynamic()) {
70 <      seleMan1_.setSelectionSet(evaluator1_.evaluate());
71 <    }else {
72 <      sprintf( painCave.errMsg,
73 <               "--sele1 must be static selection\n");
74 <      painCave.severity = OOPSE_ERROR;
75 <      painCave.isFatal = 1;
76 <      simError();
67 >    evaluator_.loadScriptString(sele);
68 >    if (!evaluator_.isDynamic()) {
69 >      seleMan_.setSelectionSet(evaluator_.evaluate());
70      }
71  
72 <    /* Set up cutoff radius and type of order parameter we are calcuating*/
72 >    // Set up cutoff radius and order of the Legendre Polynomial:
73 >
74      lNumber_ = lNumber;
75      rCut_ = rCut;
76 <    mSize_ = 2*lNumber_+1;
76 >    mSize_ = 2*lNumber_+1;    
77  
78 <    int i;
85 <    int j;
86 <    StuntDouble* sd1;
87 <    StuntDouble* sd2;
88 <    for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan1_.beginSelected(j);
89 <         sd1 != NULL && sd2 != NULL;
90 <         sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) {
91 <      for (sd2 = seleMan1_.beginSelected(j),sd2
92 <             sdPairs_.push_back(std::make_pair(sd1, sd2));
93 <           }
78 >    // Set the l for the spherical harmonic, it doesn't change
79  
80 +    sphericalHarmonic.setL(lNumber_);
81  
82 <    }
82 >    delta_Q = 1.0 / nbins;
83 >    delta_W = 2.0 / nbins;
84  
85 <    void BondOrderParameter::process
86 <      () {
100 <      Molecule* mol;
101 <      RigidBody* rb;
102 <      SimInfo::MoleculeIterator mi;
103 <      Molecule::RigidBodyIterator rbIter;
104 <      RealType theta;
105 <      RealType phi;
106 <      RealType r;
107 <      RealType dist;
108 <      RealType* QBar_lm;
109 <      RealType QSq_l;
110 <      int nBonds;
111 <      int m, m_index;
112 <      RealSphericalHarmonic sphericalHarmonic;
85 >    Q_histogram_.resize(nbins);
86 >    W_histogram_.resize(nbins);
87  
88 +  }
89  
90 <      DumpReader reader(info_, dumpFilename_);
91 <      int nFrames = reader.getNFrames();
90 >  void BondOrderParameter::initalizeHistogram() {
91 >    std::fill(Q_histogram_.begin(), Q_histogram_.end(), 0);
92 >    std::fill(W_histogram_.begin(), W_histogram_.end(), 0);
93 >  }
94  
95 <      /*Set the l for the spherical harmonic, it doesn't change*/
96 <      sphericalHarmonic.setL(lNumber_);
95 >  void BondOrderParameter::process() {
96 >    Molecule* mol;
97 >    Atom* atom;
98 >    RigidBody* rb;
99 >    SimInfo::MoleculeIterator mi;
100 >    Molecule::RigidBodyIterator rbIter;
101 >    Molecule::AtomIterator ai;
102 >    StuntDouble* sd;
103 >    RealType theta;
104 >    RealType phi;
105 >    RealType r;
106 >    RealType dist;
107 >    std::map<int, RealType> QBar_lm;
108 >    RealType QSq_l;
109 >    RealType Q_l;
110 >    int nBonds;
111 >    RealSphericalHarmonic sphericalHarmonic;
112 >    int i, j;
113 >  
114 >  
115 >    DumpReader reader(info_, dumpFilename_);    
116 >    int nFrames = reader.getNFrames();
117  
121 –      for (int i = 0; i < nFrames; i += step_) {
122 –        reader.readFrame(i);
123 –        currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
124 –        nBonds = 0;
118  
119 <        for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
120 <          //change the positions of atoms which belong to the rigidbodies
121 <          for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
122 <            rb->updateAtoms();
123 <          }
119 >    for (int istep = 0; istep < nFrames; istep += step_) {
120 >      reader.readFrame(istep);
121 >      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
122 >      
123 >      if (evaluator_.isDynamic()) {
124 >        seleMan_.setSelectionSet(evaluator_.evaluate());
125 >      }
126  
127 <        }
127 >      // update the positions of atoms which belong to the rigidbodies
128  
129 +      for (mol = info_->beginMolecule(mi); mol != NULL;
130 +           mol = info_->nextMolecule(mi)) {
131 +        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
132 +             rb = mol->nextRigidBody(rbIter)) {
133 +          rb->updateAtoms();
134 +        }        
135 +      }      
136 +      
137 +      // outer loop is over the selected StuntDoubles:
138  
139 <        /* Setup QBar */
140 <        QBar_lm = new double[mSize_];
139 >      for (sd = seleMan_.beginSelected(i); sd != NULL;
140 >           sd = seleMan_.nextSelected(i)) {
141  
142 <        /* Calculate "bonds" and build Q_lm(r) where Q_lm = Y_lm(theta(r),phi(r)) */
139 <        for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
140 <          Vector3d vec = j->first->getPos() - j->second->getPos();
141 <          currentSnapshot_->wrapVector(vec);
142 <          /* The spherical harmonics are wrt any arbitray coordiate sysetm,
143 <           * we choose standard spherical coordinates */
144 <          r = sqrt(pow(vec.x(),2)+pow(vec.y(),2)+pow(vec.z(),2));
142 >        // For this central atom, zero out nBonds and QBar_lm
143  
144 <          /* Check to see if neighbor is in bond cuttoff*/
145 <          if (r<rCut_){
146 <            theta = atan(vec.y()/vec.x());
147 <            phi = acos(vec.z()/r);
148 <            for(int m_index = 0; m_index < mSize_; m_index++){
149 <              sphericalHarmonic.setM(m_index-lNumber_);
150 <              QBar_lm(m_index) += sphericalHarmonic.getValueAt(theta,phi);
151 <            }
152 <            nBonds++;
153 <          }
154 <        }
144 >        nBonds = 0;
145 >      
146 >        for (int m = -lNumber_; m <= lNumber_; m++) {
147 >          QBar_lm[m] = 0.0;
148 >        }
149 >        
150 >        // inner loop is over all other atoms in the system:
151 >        
152 >        for (mol = info_->beginMolecule(mi); mol != NULL;
153 >             mol = info_->nextMolecule(mi)) {
154 >          for (atom = mol->beginAtom(ai); atom != NULL;
155 >               atom = mol->nextAtom(ai)) {
156  
158 –        /*Normalize Qbar by number of Bonds*/
159 –        for ( int m_index = 0;m_index < mSize_; m_index++){
160 –          QBar_lm(m_index) = QBar_lm(m_index)/nBonds;
161 –        }
157  
158 <
159 <      }
160 <
161 <      /*Normalize by number of frames*/
162 <      for ( int m_index = 0;m_index < mSize_; m_index++){
163 <        QBar_lm(m_index) = QBar_lm(m_index)/nFrames;
164 <      }
158 >            Vector3d vec = sd->getPos() - atom->getPos();      
159 >            currentSnapshot_->wrapVector(vec);
160 >            
161 >            // Calculate "bonds" and build Q_lm(r) where
162 >            //      Q_lm = Y_lm(theta(r),phi(r))                
163 >            // The spherical harmonics are wrt any arbitrary coordinate
164 >            // system, we choose standard spherical coordinates
165 >            
166 >            r = sqrt(pow(vec.x(),2)+pow(vec.y(),2)+pow(vec.z(),2));
167 >            
168 >            // Check to see if neighbor is in bond cutoff
169 >            
170 >            if (r < rCut_) {            
171 >              theta = atan2(vec.y(), vec.x());
172 >              phi = acos(vec.z()/r);
173 >              for(int m = -lNumber_; m <= lNumber_; m++){
174 >                sphericalHarmonic.setM(m);
175 >                QBar_lm[m] += sphericalHarmonic.getValueAt(theta,phi);
176 >              }
177 >              nBonds++;
178 >            }  
179 >          }
180 >        }
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 <      /* Find second order invariant Q_l*/
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_];
216 >          }
217 >        }
218  
219 <      for (int m_index = 0 ;m_index <= sizeM_; m++){
176 <        QSq_l += pow(QBar_lm(m),2);
177 <      }
178 <      Q_l_ = sqrt((4*NumericConstant::PI/lNumber_+1)*QSq_l);
219 >        W_l_hat = W_l_ / pow(QSq_l, 1.5);
220  
221 <      /* Find Third Order Invariant W_l*/
221 >        // accumulate histogram data for Q_l and W_l_hat:
222  
223 <      /* Make arrays for Wigner3jm */
224 <      double* THRCOF = new double[mSize_];
184 <      /* Variables for Wigner routine */
185 <      double l_ = (double)lNumber_;
186 <      double m2Min;
187 <      double m2Max;
188 <      int error;
189 <      int m1;
190 <      int m2;
191 <      int m3;
192 <
193 <      for (int m1 = -lNumber_;m <= lNumber_;m1++){
194 <        /* Zero work array */
195 <        for (i=0; i<mSize_;i++){
196 <          THRCOF[i] = 0.0;      
197 <        }
198 <        /* Get wigner coefficients */
199 <        Wigner3jm(&l_,&l_,&l_,&(double)m1,&m2Min,&m2Max&,THRCOF,&mSize_,&error);
200 <        for (m_index=1; i<(m2Max-M2Min-1.0);m_index++){
201 <          m2 = floor(m2Min) + m_index - 1;
202 <          m3 = -m1-m2;
203 <          W_l_ += THRCOF(m_index)*QBar_lm(m1+lNumber_)*QBar_lm(m2+lNumber_)*QBar_lm(m3+lNumber_);
204 <        }
223 >        collectHistogram(Q_l, W_l_hat);
224 >                
225        }
226 +    }
227 +    
228 +    // Normalize by number of frames
229 +    for (int m = -lNumber_; m <= lNumber_; m++){
230 +      QBar_lm[m] /=  nFrames;  
231 +    }
232 +    
233 +    
234 +    
235 +    
236 +    
237 +    writeOrderParameter();
238 +    
239 +  }
240  
241 +  
242 +  void BondOrderParameter::processHistogram() {
243 +    
244 +    int nPairs = getNPairs();
245 +    RealType volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
246 +    RealType pairDensity = nPairs /volume * 2.0;
247 +    RealType pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0;
248  
249 <      writeOrderParameter();
249 >    for(int i = 0 ; i < histogram_.size(); ++i){
250  
251 +      RealType rLower = i * deltaR_;
252 +      RealType rUpper = rLower + deltaR_;
253 +      RealType volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower );
254 +      RealType nIdeal = volSlice * pairConstant;
255 +
256 +      avgGofr_[i] += histogram_[i] / nIdeal;    
257      }
258  
259 +  }
260  
261 <    void BondOrderParameter::writeOrderParameter() {
261 >  void BondOrderParameter::collectHistogram(RealType Q_l, RealType W_l_hat) {
262  
263 <      std::ofstream os(getOutputFileName().c_str());
264 <      os << "#radial distribution function\n";
265 <      os<< "#selection1: (" << selectionScript1_ << ")\t";
266 <      os << "selection2: (" << selectionScript2_ << ")\n";
267 <      os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\n";
263 >    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 >  
273  
274 <      for (std::size_t i = 0; i < orderParams_.size(); ++i) {
222 <        os <<  orderParams_[i].p2 << "\t"
223 <           <<  orderParams_[i].director[0] << "\t"
224 <           <<  orderParams_[i].director[1] << "\t"
225 <           <<  orderParams_[i].director[2] << "\t"
226 <           <<  orderParams_[i].angle << "\n";
274 >  void BondOrderParameter::writeOrderParameter() {
275  
276 <      }
276 >    std::ofstream os(getOutputFileName().c_str());
277 >    os << "#Bond Order Parameter\n";
278 >    os << "#selection: (" << selectionScript_ << ")\n";
279  
280 <    }
280 >    for (std::size_t i = 0; i < orderParams_.size(); ++i) {
281 >      os <<  orderParams_[i].p2 << "\t"
282 >         <<  orderParams_[i].director[0] << "\t"
283 >         <<  orderParams_[i].director[1] << "\t"
284 >         <<  orderParams_[i].director[2] << "\t"
285 >         <<  orderParams_[i].angle << "\n";
286  
287 +    }
288    }
289  
290 +
291 +
292 + }

Diff Legend

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