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

Comparing trunk/src/applications/staticProps/RippleOP.cpp (file contents):
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

# Line 48 | Line 48 | namespace OpenMD {
48   #include "types/MultipoleAdapter.hpp"
49   namespace OpenMD {
50  
51 <
52 <  RippleOP::RippleOP(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2)
53 <    : StaticAnalyser(info, filename),
54 <      selectionScript1_(sele1), selectionScript2_(sele2), evaluator1_(info), evaluator2_(info),
55 <      seleMan1_(info), seleMan2_(info){
51 >  
52 >  RippleOP::RippleOP(SimInfo* info, const std::string& filename,
53 >                     const std::string& sele1, const std::string& sele2)
54 >    : StaticAnalyser(info, filename),
55 >      selectionScript1_(sele1), selectionScript2_(sele2),
56 >      seleMan1_(info), seleMan2_(info), evaluator1_(info), evaluator2_(info) {
57      
58      setOutputName(getPrefix(filename) + ".rp2");
59      
# Line 124 | Line 125 | namespace OpenMD {
125        int nTail=0;
126        RealType sumZ = 0.0;
127        
128 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
128 >      for (mol = info_->beginMolecule(mi); mol != NULL;
129 >           mol = info_->nextMolecule(mi)) {
130          //change the positions of atoms which belong to the rigidbodies
131 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
131 >        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
132 >             rb = mol->nextRigidBody(rbIter)) {
133            rb->updateAtoms();
134          }
135        }      
136  
137 <      for (sd3 = seleMan2_.beginSelected(i1); sd3 != NULL; sd3 = seleMan2_.nextSelected(i1)) {
137 >      for (sd3 = seleMan2_.beginSelected(i1); sd3 != NULL;
138 >           sd3 = seleMan2_.nextSelected(i1)) {
139          Vector3d pos1 = sd3->getPos();
140          if (usePeriodicBoundaryConditions_)
141            currentSnapshot_->wrapVector(pos1);
142          sd3->setPos(pos1);
143        }
144  
145 <      for (sd3 = seleMan2_.beginSelected(i1); sd3 != NULL; sd3 = seleMan2_.nextSelected(i1)) {
145 >      for (sd3 = seleMan2_.beginSelected(i1); sd3 != NULL;
146 >           sd3 = seleMan2_.nextSelected(i1)) {
147          Vector3d pos1 = sd3->getPos();
148          sumZ += pos1.z();
149        }
150        RealType avgZ = sumZ / (RealType) nMolecules;
151        
152 <      Mat3x3d orderTensorHeadUpper(0.0), orderTensorTail(0.0), orderTensorHeadLower(0.0);
153 <      //      for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
154 <      for (j1 = seleMan1_.beginSelected(i1); j1 != NULL; j1 = seleMan1_.nextSelected(i1)) {
152 >      Mat3x3d orderTensorHeadUpper(0.0);
153 >      Mat3x3d orderTensorTail(0.0);
154 >      Mat3x3d orderTensorHeadLower(0.0);
155 >      for (j1 = seleMan1_.beginSelected(i1); j1 != NULL;
156 >           j1 = seleMan1_.nextSelected(i1)) {
157          Vector3d pos = j1->getPos();
158          if (usePeriodicBoundaryConditions_)
159            currentSnapshot_->wrapVector(pos);
# Line 173 | Line 180 | namespace OpenMD {
180          orderTensorHeadUpper +=outProduct(vecHeadUpper, vecHeadUpper);
181          orderTensorHeadLower +=outProduct(vecHeadLower, vecHeadLower);
182        }
183 <      for (j2 = seleMan2_.beginSelected(i1); j2 != NULL; j2 = seleMan2_.nextSelected(i1)) {
183 >      for (j2 = seleMan2_.beginSelected(i1); j2 != NULL;
184 >           j2 = seleMan2_.nextSelected(i1)) {
185          // The lab frame vector corresponding to the body-fixed
186          // z-axis is simply the second column of A.transpose()
187          // or, identically, the second row of A itself.
# Line 192 | Line 200 | namespace OpenMD {
200        orderTensorTail -= (RealType)(1.0/3.0) * Mat3x3d::identity();  
201        
202        Vector3d eigenvaluesHeadUpper, eigenvaluesHeadLower, eigenvaluesTail;
203 <      Mat3x3d eigenvectorsHeadUpper, eigenvectorsHeadLower, eigenvectorsTail;    
204 <      Mat3x3d::diagonalize(orderTensorHeadUpper, eigenvaluesHeadUpper, eigenvectorsHeadUpper);
205 <      Mat3x3d::diagonalize(orderTensorHeadLower, eigenvaluesHeadLower, eigenvectorsHeadLower);
203 >      Mat3x3d eigenvectorsHeadUpper, eigenvectorsHeadLower, eigenvectorsTail;
204 >      Mat3x3d::diagonalize(orderTensorHeadUpper, eigenvaluesHeadUpper,
205 >                           eigenvectorsHeadUpper);
206 >      Mat3x3d::diagonalize(orderTensorHeadLower, eigenvaluesHeadLower,
207 >                           eigenvectorsHeadLower);
208        Mat3x3d::diagonalize(orderTensorTail, eigenvaluesTail, eigenvectorsTail);
209        
210 <      int whichUpper, whichLower, whichTail;
210 >      int whichUpper(-1), whichLower(-1), whichTail(-1);
211        RealType maxEvalUpper = 0.0;
212        RealType maxEvalLower = 0.0;
213        RealType maxEvalTail = 0.0;
# Line 223 | Line 233 | namespace OpenMD {
233        }
234        RealType p2Tail = 1.5 * maxEvalTail;
235        
236 <      //the eigen vector is already normalized in SquareMatrix3::diagonalize
236 >      //the eigenvector is already normalized in SquareMatrix3::diagonalize
237        Vector3d directorHeadUpper = eigenvectorsHeadUpper.getColumn(whichUpper);
238        if (directorHeadUpper[0] < 0) {
239          directorHeadUpper.negate();

Diff Legend

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