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 995 by chuckv, Tue Jun 27 16:36:25 2006 UTC vs.
Revision 1001 by chuckv, Fri Jul 7 14:18:36 2006 UTC

# Line 39 | Line 39
39   * such damages.
40   */
41  
42 +
43 + /* Creates orientational bond order parameters as outlined by
44 + *     Bond-orientaional order in liquids and glasses, Steinhart,Nelson,Ronchetti
45 + *     Phys Rev B, 28,784,1983
46 + *
47 + */
48 +
49   #include "applications/staticProps/BondOrderParameter.hpp"
50   #include "utils/simError.h"
51   #include "io/DumpReader.hpp"
52   #include "primitives/Molecule.hpp"
53   #include "utils/NumericConstant.hpp"
54 + #include "math/RealSphericalHarmonic.hpp"
55   namespace oopse {
56  
57  
# Line 53 | Line 61 | BondOrderParameter::BondOrderParameter(SimInfo* info,
61      selectionScript1_(sele1), evaluator1_(info),
62      seleMan1_(info){
63  
64 <    setOutputName(getPrefix(filename) + ".p2");
64 >    setOutputName(getPrefix(filename) + ".obo");
65          
66      evaluator1_.loadScriptString(sele1);
67      evaluator2_.loadScriptString(sele2);
# Line 68 | Line 76 | BondOrderParameter::BondOrderParameter(SimInfo* info,
76          simError();  
77      }
78  
79 + /* Set up cutoff radius and type of order parameter we are calcuating*/
80 +        lNumber_ = lNumber;
81 +        rCut_ = rCut;
82 +        mSize_ = 2*lNumber_+1;
83  
84    int i;
85    int j;
# Line 88 | Line 100 | void BondOrderParameter::process() {
100    RigidBody* rb;
101    SimInfo::MoleculeIterator mi;
102    Molecule::RigidBodyIterator rbIter;
103 +  RealType theta;
104 +  RealType phi;
105 +  RealType r;
106 +  RealType dist;
107 +  RealType* QBar_lm;
108 +  int nBonds;
109 +  RealSphericalHarmonic sphericalHarmonic;
110    
111 +  
112    DumpReader reader(info_, dumpFilename_);    
113    int nFrames = reader.getNFrames();
114  
115 +   /*Set the l for the spherical harmonic, it doesn't change*/
116 +        sphericalHarmonic.setL(lNumber_);
117 +
118    for (int i = 0; i < nFrames; i += step_) {
119      reader.readFrame(i);
120      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
121 <
121 >        nBonds = 0;
122      
123      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
124          //change the positions of atoms which belong to the rigidbodies
# Line 105 | Line 128 | void BondOrderParameter::process() {
128          
129      }      
130  
131 <      Mat3x3d orderTensor(0.0);
131 >      /* Calculate "bonds" and build Q_lm(r) where Q_lm = Y_lm(theta(r),phi(r)) */
132        for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
133            Vector3d vec = j->first->getPos() - j->second->getPos();
134            currentSnapshot_->wrapVector(vec);
135 <          vec.normalize();
136 <          orderTensor +=outProduct(vec, vec);
135 > /* The spherical harmonics are wrt any arbitray coordiate sysetm,
136 >  * we choose standard spherical coordinates */
137 >                r = sqrt(pow(vec.x(),2)+pow(vec.y(),2)+pow(vec.z(),2));
138 >                
139 >                /* Check to see if neighbor is in bond cuttoff*/
140 >                if (r<rCut_){          
141 >            theta = atan(vec.y()/vec.x());
142 >            phi = acos(vec.z()/r);
143 >            for(int m = -lNumber_; m <= lNumber_; m++){
144 >                sphericalHarmonic.setM(m);
145 >                QBar_lm(m) += sphericalHarmonic.getValueAt(theta,phi);
146 >            }
147 >            nBonds++;
148 >                }  
149        }
150        
151 <      orderTensor /= sdPairs_.size();
152 <      orderTensor -= (RealType)(1.0/3.0) * Mat3x3d::identity();  
151 >     /*Normalize Qbar*/
152 >     for (int m = -lNumber_;m <= lNumber_; m++){
153 >      QBar_lm(m) = QBar_lm(m)/nBonds;  
154 >     }
155        
156 <      Vector3d eigenvalues;
157 <      Mat3x3d eigenvectors;    
158 <      Mat3x3d::diagonalize(orderTensor, eigenvalues, eigenvectors);
159 <      
160 <      int which;
161 <      RealType maxEval = 0.0;
162 <      for(int k = 0; k< 3; k++){
163 <        if(fabs(eigenvalues[k]) > maxEval){
164 <          which = k;
165 <          maxEval = fabs(eigenvalues[k]);
166 <        }
167 <      }
168 <      RealType p2 = 1.5 * maxEval;
169 <      
170 <      //the eigen vector is already normalized in SquareMatrix3::diagonalize
171 <      Vector3d director = eigenvectors.getColumn(which);
172 <      if (director[0] < 0) {
173 <          director.negate();
174 <      }  
156 >    
157 >  }
158 >  
159 >  /*Normalize by number of frames*/
160 >    for (int m = -lNumber_;m <= lNumber_; m++){
161 >      QBar_lm(m) = QBar_lm(m)/nFrames;  
162 >     }
163 >    
164 >    
165 >    
166 >     /* Find second order invariant Q_l*/
167 >    
168 >       for (int m = -lNumber_;m <= lNumber_; m++){
169 >         QSq_l += pow(QBar_lm(m),2);    
170 >        }
171 >     Q_l = sqrt((4*NumericConstant::PI/lNumber_+1)*QSq_l);
172 >    
173 >     /* Find Third Order Invariant W_l*/
174 >     for (int m = -lNumber_;m<= lNumber_;m++){
175 >        
176 >        
177 >     }
178 >    
179 >    
180 >  writeOrderParameter();
181 >  
182 > }
183  
184 <      RealType angle = 0.0;
185 <      for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
186 <          Vector3d vec = j->first->getPos() - j->second->getPos();
142 <          currentSnapshot_->wrapVector(vec);
143 <          vec.normalize();
184 > void BondOrderParameter::initalizeHistogram() {
185 >    std::fill(histogram_.begin(), histogram_.end(), 0);
186 >  }
187  
188 <          angle += acos(dot(vec, director)) ;
146 <      }
147 <      angle = angle / (sdPairs_.size() * NumericConstant::PI) * 180.0;
188 > void BondOrderParameter::processHistogram() {
189  
190 <       OrderParam param;
191 <       param.p2 = p2;
192 <       param.director = director;
193 <       param.angle = angle;
190 >    int nPairs = getNPairs();
191 >    RealType volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
192 >    RealType pairDensity = nPairs /volume * 2.0;
193 >    RealType pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0;
194  
195 <        orderParams_.push_back(param);      
195 >    for(int i = 0 ; i < histogram_.size(); ++i){
196 >
197 >      RealType rLower = i * deltaR_;
198 >      RealType rUpper = rLower + deltaR_;
199 >      RealType volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower );
200 >      RealType nIdeal = volSlice * pairConstant;
201 >
202 >      avgGofr_[i] += histogram_[i] / nIdeal;    
203 >    }
204 >
205 >  }
206 >
207 >  void BondOrderParameter::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
208 >
209 >    if (sd1 == sd2) {
210 >      return;
211 >    }
212      
213 +    Vector3d pos1 = sd1->getPos();
214 +    Vector3d pos2 = sd2->getPos();
215 +    Vector3d r12 = pos2 - pos1;
216 +    currentSnapshot_->wrapVector(r12);
217 +
218 +    RealType distance = r12.length();
219 +
220 +    if (distance < len_) {
221 +      int whichBin = distance / deltaR_;
222 +      histogram_[whichBin] += 2;
223 +    }
224    }
225  
158 –  writeP2();
159 –  
160 – }
226  
227 +
228 +
229 +
230 +
231   void BondOrderParameter::writeOrderParameter() {
232  
233      std::ofstream os(getOutputFileName().c_str());

Diff Legend

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