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

Comparing trunk/src/applications/staticProps/P2OrderParameter.cpp (file contents):
Revision 543 by tim, Thu May 26 22:45:00 2005 UTC vs.
Revision 963 by tim, Wed May 17 21:51:42 2006 UTC

# Line 43 | Line 43
43   #include "utils/simError.h"
44   #include "io/DumpReader.hpp"
45   #include "primitives/Molecule.hpp"
46 <
46 > #include "utils/NumericConstant.hpp"
47   namespace oopse {
48  
49  
# Line 112 | Line 112 | void P2OrderParameter::process() {
112    for (int i = 0; i < nFrames; i += step_) {
113      reader.readFrame(i);
114      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
115 <    Mat3x3d orderTensor(0.0);
115 >
116      
117      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
118          //change the positions of atoms which belong to the rigidbodies
# Line 122 | Line 122 | void P2OrderParameter::process() {
122          
123      }      
124  
125 +      Mat3x3d orderTensor(0.0);
126        for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
127            Vector3d vec = j->first->getPos() - j->second->getPos();
128 +          currentSnapshot_->wrapVector(vec);
129            vec.normalize();
130            orderTensor +=outProduct(vec, vec);
131        }
132        
133        orderTensor /= sdPairs_.size();
134 <      orderTensor -= 1.0/3.0 * Mat3x3d::identity();  
134 >      orderTensor -= (RealType)(1.0/3.0) * Mat3x3d::identity();  
135        
136        Vector3d eigenvalues;
137        Mat3x3d eigenvectors;    
138        Mat3x3d::diagonalize(orderTensor, eigenvalues, eigenvectors);
139        
140        int which;
141 <      double maxEval = 0.0;
141 >      RealType maxEval = 0.0;
142        for(int k = 0; k< 3; k++){
143          if(fabs(eigenvalues[k]) > maxEval){
144            which = k;
145            maxEval = fabs(eigenvalues[k]);
146          }
147        }
148 <      double p2 = 1.5 * maxEval;
148 >      RealType p2 = 1.5 * maxEval;
149        
150        //the eigen vector is already normalized in SquareMatrix3::diagonalize
151        Vector3d director = eigenvectors.getColumn(which);
# Line 151 | Line 153 | void P2OrderParameter::process() {
153            director.negate();
154        }  
155  
156 <      double angle = 0.0;
156 >      RealType angle = 0.0;
157        for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
158            Vector3d vec = j->first->getPos() - j->second->getPos();
159 +          currentSnapshot_->wrapVector(vec);
160            vec.normalize();
161  
162            angle += acos(dot(vec, director)) ;
163        }
164 <      angle /= sdPairs_.size();
164 >      angle = angle / (sdPairs_.size() * NumericConstant::PI) * 180.0;
165  
166         OrderParam param;
167         param.p2 = p2;
# Line 169 | Line 172 | void P2OrderParameter::process() {
172      
173    }
174  
175 <  writeOrderParam();
175 >  writeP2();
176    
177   }
178  
179 < void P2OrderParameter::writeOrderParam() {
179 > void P2OrderParameter::writeP2() {
180  
181      std::ofstream os(getOutputFileName().c_str());
182      os << "#radial distribution function\n";

Diff Legend

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