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 1817 by gezelter, Wed Dec 12 19:05:43 2012 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC

# Line 35 | Line 35
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39   * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40   * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
# Line 51 | Line 51 | namespace OpenMD {
51  
52    P2OrderParameter::P2OrderParameter(SimInfo* info, const string& filename,
53                                       const string& sele1)
54 <    : StaticAnalyser(info, filename), doVect_(true),
55 <      selectionScript1_(sele1), evaluator1_(info),
56 <      evaluator2_(info), seleMan1_(info), seleMan2_(info) {
54 >  : StaticAnalyser(info, filename), doVect_(true), doOffset_(false),
55 >    selectionScript1_(sele1), evaluator1_(info),
56 >    evaluator2_(info), seleMan1_(info), seleMan2_(info) {
57      
58      setOutputName(getPrefix(filename) + ".p2");
59      
# Line 62 | Line 62 | namespace OpenMD {
62  
63    P2OrderParameter::P2OrderParameter(SimInfo* info, const string& filename,
64                                       const string& sele1, const string& sele2)
65 <    : StaticAnalyser(info, filename), doVect_(false),
66 <      selectionScript1_(sele1), selectionScript2_(sele2), evaluator1_(info),
67 <      evaluator2_(info), seleMan1_(info), seleMan2_(info) {
65 >  : StaticAnalyser(info, filename), doVect_(false), doOffset_(false),
66 >    selectionScript1_(sele1), selectionScript2_(sele2), evaluator1_(info),
67 >    evaluator2_(info), seleMan1_(info), seleMan2_(info) {
68      
69      setOutputName(getPrefix(filename) + ".p2");
70      
# Line 72 | Line 72 | namespace OpenMD {
72      evaluator2_.loadScriptString(sele2);    
73    }
74  
75 +  P2OrderParameter::P2OrderParameter(SimInfo* info, const string& filename,
76 +                                     const string& sele1, int seleOffset)
77 +  : StaticAnalyser(info, filename), doVect_(false), doOffset_(true),
78 +    seleOffset_(seleOffset), selectionScript1_(sele1), evaluator1_(info),
79 +    evaluator2_(info), seleMan1_(info), seleMan2_(info) {
80 +    
81 +    setOutputName(getPrefix(filename) + ".p2");
82 +    
83 +    evaluator1_.loadScriptString(sele1);
84 +  }
85 +
86    void P2OrderParameter::process() {
87      Molecule* mol;
88      RigidBody* rb;
# Line 102 | Line 113 | namespace OpenMD {
113        Mat3x3d orderTensor(0.0);
114        vecCount = 0;
115  
116 <      if  (evaluator1_.isDynamic())
106 <        seleMan1_.setSelectionSet(evaluator1_.evaluate());
116 >      seleMan1_.setSelectionSet(evaluator1_.evaluate());
117        
118        if (doVect_) {
119          
# Line 121 | Line 131 | namespace OpenMD {
131  
132        } else {
133  
134 <        if (evaluator2_.isDynamic())
125 <          seleMan2_.setSelectionSet(evaluator2_.evaluate());
126 <            
127 <        if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) {
128 <          sprintf( painCave.errMsg,
129 <                   "In frame %d, the number of selected StuntDoubles are not\n"
130 <                   "\tthe same in --sele1 and sele2\n", i);
131 <          painCave.severity = OPENMD_INFO;
132 <          painCave.isFatal = 0;
133 <          simError();            
134 <        }
135 <        
136 <        for (sd1 = seleMan1_.beginSelected(ii),
137 <               sd2 = seleMan2_.beginSelected(jj);
138 <             sd1 != NULL && sd2 != NULL;
139 <             sd1 = seleMan1_.nextSelected(ii),
140 <               sd2 = seleMan2_.nextSelected(jj)) {
134 >        if (doOffset_) {
135  
136 <          Vector3d vec = sd1->getPos() - sd2->getPos();
136 >          for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
137 >               sd1 = seleMan1_.nextSelected(ii)) {
138  
139 <          if (usePeriodicBoundaryConditions_)
140 <            currentSnapshot_->wrapVector(vec);
139 >            // This will require careful rewriting if StaticProps is
140 >            // ever parallelized.  For an example, see
141 >            // Thermo::getTaggedAtomPairDistance
142  
143 <          vec.normalize();
143 >            int sd2Index = sd1->getGlobalIndex() + seleOffset_;
144 >            sd2 = info_->getIOIndexToIntegrableObject(sd2Index);
145 >            
146 >            Vector3d vec = sd1->getPos() - sd2->getPos();
147 >            
148 >            if (usePeriodicBoundaryConditions_)
149 >              currentSnapshot_->wrapVector(vec);
150 >            
151 >            vec.normalize();
152 >            
153 >            orderTensor +=outProduct(vec, vec);
154 >            vecCount++;
155 >          }
156 >          
157 >          orderTensor /= vecCount;
158 >        } else {
159 >          
160 >          seleMan2_.setSelectionSet(evaluator2_.evaluate());
161 >          
162 >          if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) {
163 >            sprintf( painCave.errMsg,
164 >                     "In frame %d, the number of selected StuntDoubles are\n"
165 >                     "\tnot the same in --sele1 and sele2\n", i);
166 >            painCave.severity = OPENMD_INFO;
167 >            painCave.isFatal = 0;
168 >            simError();            
169 >          }
170 >          
171 >          for (sd1 = seleMan1_.beginSelected(ii),
172 >                 sd2 = seleMan2_.beginSelected(jj);
173 >               sd1 != NULL && sd2 != NULL;
174 >               sd1 = seleMan1_.nextSelected(ii),
175 >                 sd2 = seleMan2_.nextSelected(jj)) {
176 >            
177 >            Vector3d vec = sd1->getPos() - sd2->getPos();
178 >            
179 >            if (usePeriodicBoundaryConditions_)
180 >              currentSnapshot_->wrapVector(vec);
181  
182 <          orderTensor +=outProduct(vec, vec);
183 <          vecCount++;
182 >            vec.normalize();
183 >            
184 >            orderTensor +=outProduct(vec, vec);
185 >            vecCount++;
186 >          }
187 >          
188 >          orderTensor /= vecCount;
189          }
152 –      
153 –        orderTensor /= vecCount;
190        }
191        
192        if (vecCount == 0) {
# Line 201 | Line 237 | namespace OpenMD {
237          angle = angle/(vecCount*NumericConstant::PI)*180.0;
238          
239        } else {
240 <        for (sd1 = seleMan1_.beginSelected(ii),
241 <               sd2 = seleMan2_.beginSelected(jj);
242 <             sd1 != NULL && sd2 != NULL;
243 <             sd1 = seleMan1_.nextSelected(ii),
244 <               sd2 = seleMan2_.nextSelected(jj)) {
245 <          
246 <          Vector3d vec = sd1->getPos() - sd2->getPos();
247 <          if (usePeriodicBoundaryConditions_)
248 <            currentSnapshot_->wrapVector(vec);
249 <          vec.normalize();          
250 <          angle += acos(dot(vec, director)) ;
251 <          vecCount++;
240 >        if (doOffset_) {
241 >
242 >          for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
243 >               sd1 = seleMan1_.nextSelected(ii)) {
244 >            
245 >            // This will require careful rewriting if StaticProps is
246 >            // ever parallelized.  For an example, see
247 >            // Thermo::getTaggedAtomPairDistance
248 >            
249 >            int sd2Index = sd1->getGlobalIndex() + seleOffset_;
250 >            sd2 = info_->getIOIndexToIntegrableObject(sd2Index);
251 >            
252 >            Vector3d vec = sd1->getPos() - sd2->getPos();
253 >            if (usePeriodicBoundaryConditions_)
254 >              currentSnapshot_->wrapVector(vec);
255 >            vec.normalize();          
256 >            angle += acos(dot(vec, director)) ;
257 >            vecCount++;
258 >          }
259 >          angle = angle / (vecCount * NumericConstant::PI) * 180.0;
260 >
261 >        } else {
262 >
263 >          for (sd1 = seleMan1_.beginSelected(ii),
264 >                 sd2 = seleMan2_.beginSelected(jj);
265 >               sd1 != NULL && sd2 != NULL;
266 >               sd1 = seleMan1_.nextSelected(ii),
267 >                 sd2 = seleMan2_.nextSelected(jj)) {
268 >            
269 >            Vector3d vec = sd1->getPos() - sd2->getPos();
270 >            if (usePeriodicBoundaryConditions_)
271 >              currentSnapshot_->wrapVector(vec);
272 >            vec.normalize();          
273 >            angle += acos(dot(vec, director)) ;
274 >            vecCount++;
275 >          }
276 >          angle = angle / (vecCount * NumericConstant::PI) * 180.0;
277          }
217 –        angle = angle / (vecCount * NumericConstant::PI) * 180.0;
278        }
279  
280        OrderParam param;

Diff Legend

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