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 1657 by gezelter, Tue Oct 18 13:44:44 2011 UTC vs.
Revision 1819 by gezelter, Thu Dec 13 16:57:39 2012 UTC

# Line 36 | Line 36
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).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include "applications/staticProps/P2OrderParameter.hpp"
# Line 50 | 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      
60      evaluator1_.loadScriptString(sele1);
60 –    
61 –    if (!evaluator1_.isDynamic()) {
62 –      seleMan1_.setSelectionSet(evaluator1_.evaluate());
63 –    }
61    }
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      
71      evaluator1_.loadScriptString(sele1);
72 <    evaluator2_.loadScriptString(sele2);
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 <    if (!evaluator1_.isDynamic()) {
78 <      seleMan1_.setSelectionSet(evaluator1_.evaluate());
79 <    }else {
80 <      sprintf( painCave.errMsg,
81 <               "--sele1 must be static selection\n");
82 <      painCave.severity = OPENMD_ERROR;
83 <      painCave.isFatal = 1;
84 <      simError();  
85 <    }
81 >    setOutputName(getPrefix(filename) + ".p2");
82      
83 <    if (!evaluator2_.isDynamic()) {
88 <      seleMan2_.setSelectionSet(evaluator2_.evaluate());
89 <    }else {
90 <      sprintf( painCave.errMsg,
91 <               "--sele2 must be static selection\n");
92 <      painCave.severity = OPENMD_ERROR;
93 <      painCave.isFatal = 1;
94 <      simError();  
95 <    }
96 <    
97 <    if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) {
98 <      sprintf( painCave.errMsg,
99 <               "The number of selected Stuntdoubles are not the same in --sele1 and sele2\n");
100 <      painCave.severity = OPENMD_ERROR;
101 <      painCave.isFatal = 1;
102 <      simError();  
103 <      
104 <    }
105 <    
106 <    int i;
107 <    int j;
108 <    StuntDouble* sd1;
109 <    StuntDouble* sd2;
110 <    for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan2_.beginSelected(j);
111 <         sd1 != NULL && sd2 != NULL;
112 <         sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) {
113 <
114 <      sdPairs_.push_back(make_pair(sd1, sd2));
115 <    }
83 >    evaluator1_.loadScriptString(sele1);
84    }
85  
86    void P2OrderParameter::process() {
# Line 120 | Line 88 | namespace OpenMD {
88      RigidBody* rb;
89      SimInfo::MoleculeIterator mi;
90      Molecule::RigidBodyIterator rbIter;
91 <    StuntDouble* sd;
92 <    int i, ii;
93 <  
91 >    StuntDouble* sd1;
92 >    StuntDouble* sd2;
93 >    int ii;
94 >    int jj;
95 >    int vecCount;
96 >
97      DumpReader reader(info_, dumpFilename_);    
98      int nFrames = reader.getNFrames();
99  
# Line 140 | Line 111 | namespace OpenMD {
111        }      
112  
113        Mat3x3d orderTensor(0.0);
114 +      vecCount = 0;
115  
116 +      seleMan1_.setSelectionSet(evaluator1_.evaluate());
117 +      
118        if (doVect_) {
119          
120 <        if  (evaluator1_.isDynamic())
121 <          seleMan1_.setSelectionSet(evaluator1_.evaluate());
122 <        
123 <        for (sd = seleMan1_.beginSelected(ii); sd != NULL;
150 <             sd = seleMan1_.nextSelected(ii)) {
151 <          if (sd->isDirectional()) {
152 <            Vector3d vec = sd->getA().getColumn(2);
120 >        for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
121 >             sd1 = seleMan1_.nextSelected(ii)) {
122 >          if (sd1->isDirectional()) {
123 >            Vector3d vec = sd1->getA().getColumn(2);
124              vec.normalize();
125              orderTensor += outProduct(vec, vec);
126 +            vecCount++;
127            }
128          }
129    
130 <        orderTensor /= seleMan1_.getSelectionCount();
130 >        orderTensor /= vecCount;
131  
132        } else {
133 +
134 +        if (doOffset_) {
135 +
136 +          for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
137 +               sd1 = seleMan1_.nextSelected(ii)) {
138 +
139 +            // This will require careful rewriting if StaticProps is
140 +            // ever parallelized.  For an example, see
141 +            // Thermo::getTaggedAtomPairDistance
142 +
143 +            int sd2Index = sd1->getGlobalIndex() + seleOffset_;
144 +            sd2 = info_->getIOIndexToIntegrableObject(sd2Index);
145              
146 <        for (vector<pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin();
147 <             j != sdPairs_.end(); ++j) {
148 <          Vector3d vec = j->first->getPos() - j->second->getPos();
149 <          if (usePeriodicBoundaryConditions_)
150 <            currentSnapshot_->wrapVector(vec);
151 <          vec.normalize();
152 <          orderTensor +=outProduct(vec, vec);
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 >            vec.normalize();
183 >            
184 >            orderTensor +=outProduct(vec, vec);
185 >            vecCount++;
186 >          }
187 >          
188 >          orderTensor /= vecCount;
189          }
170 –      
171 –        orderTensor /= sdPairs_.size();
190        }
191        
192 +      if (vecCount == 0) {
193 +          sprintf( painCave.errMsg,
194 +                   "In frame %d, the number of selected vectors was zero.\n"
195 +                   "\tThis will not give a meaningful order parameter.", i);
196 +          painCave.severity = OPENMD_ERROR;
197 +          painCave.isFatal = 1;
198 +          simError();        
199 +      }
200  
201        orderTensor -= (RealType)(1.0/3.0) * Mat3x3d::identity();  
202        
# Line 196 | Line 222 | namespace OpenMD {
222        }  
223  
224        RealType angle = 0.0;
225 <      RealType p4 = 0.0;
225 >      vecCount = 0;
226        
227        if (doVect_) {
228 <        for (sd = seleMan1_.beginSelected(ii); sd != NULL;
229 <             sd = seleMan1_.nextSelected(ii)) {
230 <          if (sd->isDirectional()) {
231 <            Vector3d vec = sd->getA().getColumn(2);
228 >        for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
229 >             sd1 = seleMan1_.nextSelected(ii)) {
230 >          if (sd1->isDirectional()) {
231 >            Vector3d vec = sd1->getA().getColumn(2);
232              vec.normalize();
233 <            RealType myDot = dot(vec, director);
234 <            angle += acos(myDot);
209 <            p4 += (35.0 * pow(myDot, 4) - 30.0 * pow(myDot, 2) + 3.0) / 8.0;
233 >            angle += acos(dot(vec, director));
234 >            vecCount++;
235            }
236          }
237 <        angle = angle/(seleMan1_.getSelectionCount()*NumericConstant::PI)*180.0;
238 <        p4 /= (seleMan1_.getSelectionCount());
237 >        angle = angle/(vecCount*NumericConstant::PI)*180.0;
238 >        
239        } else {
240 <        for (vector<pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
241 <          Vector3d vec = j->first->getPos() - j->second->getPos();
242 <          if (usePeriodicBoundaryConditions_)
243 <            currentSnapshot_->wrapVector(vec);
244 <          vec.normalize();          
245 <          RealType myDot = dot(vec, director);
246 <          angle += acos(myDot);
247 <          p4 += (35.0 * pow(myDot, 4) - 30.0 * pow(myDot, 2) + 3.0) / 8.0;
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          }
224 –        angle = angle / (sdPairs_.size() * NumericConstant::PI) * 180.0;
225 –        p4 /= (seleMan1_.getSelectionCount());
278        }
279  
280        OrderParam param;
281        param.p2 = p2;
282        param.director = director;
283        param.angle = angle;
232 –      param.p4 = p4;
284  
285        orderParams_.push_back(param);      
286      
# Line 247 | Line 298 | namespace OpenMD {
298      if (!doVect_) {
299        os << "selection2: (" << selectionScript2_ << ")\n";
300      }
301 <    os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\t<p4>\n";    
301 >    os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\n";    
302  
303      for (size_t i = 0; i < orderParams_.size(); ++i) {
304        os <<  orderParams_[i].p2 << "\t"
305           <<  orderParams_[i].director[0] << "\t"
306           <<  orderParams_[i].director[1] << "\t"
307           <<  orderParams_[i].director[2] << "\t"
308 <         <<  orderParams_[i].angle << "\t"
258 <         <<  orderParams_[i].p4 << "\n";
308 >         <<  orderParams_[i].angle << "\n";
309  
310      }
311  

Diff Legend

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