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 1040 by gezelter, Tue Sep 19 21:20:15 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;
79 <    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 >    deltaQ_ = 1.0 / nbins;
79 >    deltaW_ = 2.0 / nbins;
80  
81 +    Q_histogram_.resize(nbins);
82 +    W_histogram_.resize(nbins);
83  
84 <    }
84 >  }
85  
86 <    void BondOrderParameter::process
87 <      () {
88 <      Molecule* mol;
89 <      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;
86 >  void BondOrderParameter::initalizeHistogram() {
87 >    std::fill(Q_histogram_.begin(), Q_histogram_.end(), 0);
88 >    std::fill(W_histogram_.begin(), W_histogram_.end(), 0);
89 >  }
90  
91 +  void BondOrderParameter::process() {
92 +    Molecule* mol;
93 +    Atom* atom;
94 +    RigidBody* rb;
95 +    SimInfo::MoleculeIterator mi;
96 +    Molecule::RigidBodyIterator rbIter;
97 +    Molecule::AtomIterator ai;
98 +    StuntDouble* sd;
99 +    RealType theta;
100 +    RealType phi;
101 +    RealType r;
102 +    RealType dist;
103 +    std::map<int, RealType> QBar_lm;
104 +    RealType QSq_l;
105 +    RealType Q_l;
106 +    int nBonds;
107 +    RealSphericalHarmonic sphericalHarmonic;
108 +    int i, j;
109 +  
110 +    // Set the l for the spherical harmonic, it doesn't change
111 +    sphericalHarmonic.setL(lNumber_);
112  
115 –      DumpReader reader(info_, dumpFilename_);
116 –      int nFrames = reader.getNFrames();
113  
114 <      /*Set the l for the spherical harmonic, it doesn't change*/
115 <      sphericalHarmonic.setL(lNumber_);
114 >    DumpReader reader(info_, dumpFilename_);    
115 >    int nFrames = reader.getNFrames();
116  
121 –      for (int i = 0; i < nFrames; i += step_) {
122 –        reader.readFrame(i);
123 –        currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
124 –        nBonds = 0;
117  
118 <        for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
119 <          //change the positions of atoms which belong to the rigidbodies
120 <          for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
121 <            rb->updateAtoms();
122 <          }
118 >    for (int istep = 0; istep < nFrames; istep += step_) {
119 >      reader.readFrame(istep);
120 >      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
121 >      
122 >      if (evaluator_.isDynamic()) {
123 >        seleMan_.setSelectionSet(evaluator_.evaluate());
124 >      }
125  
126 <        }
126 >      // update the positions of atoms which belong to the rigidbodies
127  
128 +      for (mol = info_->beginMolecule(mi); mol != NULL;
129 +           mol = info_->nextMolecule(mi)) {
130 +        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
131 +             rb = mol->nextRigidBody(rbIter)) {
132 +          rb->updateAtoms();
133 +        }        
134 +      }      
135 +      
136 +      // outer loop is over the selected StuntDoubles:
137  
138 <        /* Setup QBar */
139 <        QBar_lm = new double[mSize_];
138 >      for (sd = seleMan_.beginSelected(i); sd != NULL;
139 >           sd = seleMan_.nextSelected(i)) {
140  
141 <        /* 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));
141 >        // For this central atom, zero out nBonds and QBar_lm
142  
143 <          /* Check to see if neighbor is in bond cuttoff*/
144 <          if (r<rCut_){
145 <            theta = atan(vec.y()/vec.x());
146 <            phi = acos(vec.z()/r);
147 <            for(int m_index = 0; m_index < mSize_; m_index++){
148 <              sphericalHarmonic.setM(m_index-lNumber_);
149 <              QBar_lm(m_index) += sphericalHarmonic.getValueAt(theta,phi);
150 <            }
151 <            nBonds++;
152 <          }
153 <        }
143 >        nBonds = 0;
144 >      
145 >        for (int m = -lNumber_; m <= lNumber_; m++) {
146 >          QBar_lm[m] = 0.0;
147 >        }
148 >        
149 >        // inner loop is over all other atoms in the system:
150 >        
151 >        for (mol = info_->beginMolecule(mi); mol != NULL;
152 >             mol = info_->nextMolecule(mi)) {
153 >          for (atom = mol->beginAtom(ai); atom != NULL;
154 >               atom = mol->nextAtom(ai)) {
155  
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 –        }
156  
157 <
158 <      }
159 <
160 <      /*Normalize by number of frames*/
161 <      for ( int m_index = 0;m_index < mSize_; m_index++){
162 <        QBar_lm(m_index) = QBar_lm(m_index)/nFrames;
163 <      }
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 >        }
180 >        
181 >        // Normalize Qbar
182 >        for (int m = -lNumber_;m <= lNumber_; m++){
183 >          QBar_lm[m] /= nBonds;
184 >        }
185  
186 +        // Find second order invariant Q_l
187  
188 +        QSq_l = 0.0;
189 +        for (int m = -lNumber_; m <= lNumber_; m++){
190 +          QSq_l += pow(QBar_lm[m], 2);  
191 +        }
192 +        Q_l = sqrt(QSq_l*(4.0 * NumericConstant::PI / (2.0*(RealType)lNumber_ + 1)));
193 +    
194 +        // Find Third Order Invariant W_l
195  
196 <      /* Find second order invariant Q_l*/
196 >        // Make arrays for Wigner3jm
197 >        double* THRCOF = new double[mSize_];
198 >        // Variables for Wigner routine
199 >        double l_ = (double)lNumber_;
200 >        double m2Min, m2Max;
201 >        int error, m1, m2, m3;
202  
203 <      for (int m_index = 0 ;m_index <= sizeM_; m++){
204 <        QSq_l += pow(QBar_lm(m),2);
205 <      }
206 <      Q_l_ = sqrt((4*NumericConstant::PI/lNumber_+1)*QSq_l);
203 >        RealType W_l;
204 >        RealType W_l_hat;
205 >        W_l = 0.0;
206 >        for (int m1 = -lNumber_; m1 <= lNumber_; m1++) {
207 >          // Zero work array
208 >          for (int ii = 0; ii < mSize_; ii++){
209 >            THRCOF[i] = 0.0;
210 >          }
211 >          // Get Wigner coefficients
212 >          Wigner3jm(&l_, &l_, &l_, &(double)m1, &m2Min, &m2Max, THRCOF, &mSize_, &error);
213 >          for (int m_index = 1; i < (int)(m2Max - m2Min-1.0); m_index++) {
214 >            m2 = floor(m2Min) + m_index - 1;
215 >            m3 = -m1-m2;
216 >            W_l += THRCOF[m_index]*QBar_lm[m1+lNumber_]*QBar_lm[m2+lNumber_]*QBar_lm[m3+lNumber_];
217 >          }
218 >        }
219  
220 <      /* Find Third Order Invariant W_l*/
220 >        W_l_hat = W_l / pow(QSq_l, 1.5);
221  
222 <      /* Make arrays for Wigner3jm */
183 <      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;
222 >        // accumulate histogram data for Q_l and W_l_hat:
223  
224 <      for (int m1 = -lNumber_;m <= lNumber_;m1++){
225 <        /* 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 <        }
224 >        collectHistogram(Q_l, W_l_hat);
225 >                
226        }
227 +    }
228 +    
229 +    // Normalize by number of frames
230 +    for (int m = -lNumber_; m <= lNumber_; m++){
231 +      QBar_lm[m] /=  nFrames;  
232 +    }
233 +    
234 +    
235 +    
236 +    
237 +    
238 +    writeOrderParameter();
239 +    
240 +  }
241  
242 +  
243 +  void BondOrderParameter::processHistogram() {
244 +    
245 +    int nPairs = getNPairs();
246 +    RealType volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
247 +    RealType pairDensity = nPairs /volume * 2.0;
248 +    RealType pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0;
249  
250 <      writeOrderParameter();
250 >    for(int i = 0 ; i < histogram_.size(); ++i){
251  
252 +      RealType rLower = i * deltaR_;
253 +      RealType rUpper = rLower + deltaR_;
254 +      RealType volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower );
255 +      RealType nIdeal = volSlice * pairConstant;
256 +
257 +      avgGofr_[i] += histogram_[i] / nIdeal;    
258      }
259  
260 +  }
261  
262 <    void BondOrderParameter::writeOrderParameter() {
262 >  void BondOrderParameter::collectHistogram(RealType Q_l, RealType W_l_hat) {
263  
264 <      std::ofstream os(getOutputFileName().c_str());
265 <      os << "#radial distribution function\n";
266 <      os<< "#selection1: (" << selectionScript1_ << ")\t";
267 <      os << "selection2: (" << selectionScript2_ << ")\n";
268 <      os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\n";
264 >    if (Q_l < MaxQ_) {
265 >      int whichBin = Q_l / deltaQ_;
266 >      Q_histogram_[whichBin] += 1;
267 >    }
268 >    if (W_l_hat < MaxW_) {
269 >      int whichBin = W_l_hat / deltaW_;
270 >      W_histogram_[whichBin] += 1;
271 >    }
272 >  }
273 >  
274  
275 <      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";
275 >  void BondOrderParameter::writeOrderParameter() {
276  
277 <      }
277 >    std::ofstream os(getOutputFileName().c_str());
278 >    os << "#Bond Order Parameter\n";
279 >    os << "#selection: (" << selectionScript_ << ")\n";
280  
281 <    }
281 >    for (std::size_t i = 0; i < orderParams_.size(); ++i) {
282 >      os <<  orderParams_[i].p2 << "\t"
283 >         <<  orderParams_[i].director[0] << "\t"
284 >         <<  orderParams_[i].director[1] << "\t"
285 >         <<  orderParams_[i].director[2] << "\t"
286 >         <<  orderParams_[i].angle << "\n";
287  
288 +    }
289    }
290  
291 +
292 +
293 + }

Diff Legend

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