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

Comparing trunk/src/applications/staticProps/pAngle.cpp (file contents):
Revision 1796 by gezelter, Mon Sep 10 18:38:44 2012 UTC vs.
Revision 1991 by gezelter, Wed Apr 23 20:34:17 2014 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   * [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
41   */
# Line 53 | Line 53 | namespace OpenMD {
53   namespace OpenMD {
54    
55    pAngle::pAngle(SimInfo* info, const std::string& filename,
56 <                 const std::string& sele, int nthetabins)
57 <    : StaticAnalyser(info, filename), selectionScript_(sele),  
58 <      evaluator_(info), seleMan_(info), nThetaBins_(nthetabins){
56 >                 const std::string& sele1, int nthetabins)
57 >    : StaticAnalyser(info, filename), selectionScript1_(sele1),  
58 >      evaluator1_(info),  evaluator2_(info), seleMan1_(info), seleMan2_(info),
59 >      nThetaBins_(nthetabins),
60 >      doVect_(true), doOffset_(false) {
61      
62 <    evaluator_.loadScriptString(sele);
63 <    if (!evaluator_.isDynamic()) {
64 <      seleMan_.setSelectionSet(evaluator_.evaluate());
62 >    setOutputName(getPrefix(filename) + ".pAngle");
63 >
64 >    evaluator1_.loadScriptString(sele1);
65 >    if (!evaluator1_.isDynamic()) {
66 >      seleMan1_.setSelectionSet(evaluator1_.evaluate());
67      }            
68      
69      count_.resize(nThetaBins_);
70 <    histogram_.resize(nThetaBins_);
70 >    histogram_.resize(nThetaBins_);
71 >  }
72 >
73 >  pAngle::pAngle(SimInfo* info, const std::string& filename,
74 >                 const std::string& sele1, const std::string& sele2,
75 >                 int nthetabins)
76 >    : StaticAnalyser(info, filename), selectionScript1_(sele1),
77 >      selectionScript2_(sele2), evaluator1_(info), evaluator2_(info),
78 >      seleMan1_(info), seleMan2_(info), nThetaBins_(nthetabins),
79 >      doVect_(false), doOffset_(false) {
80      
81      setOutputName(getPrefix(filename) + ".pAngle");
82 +
83 +    evaluator1_.loadScriptString(sele1);
84 +    if (!evaluator1_.isDynamic()) {
85 +      seleMan1_.setSelectionSet(evaluator1_.evaluate());
86 +    }            
87 +    
88 +    evaluator2_.loadScriptString(sele2);
89 +    if (!evaluator2_.isDynamic()) {
90 +      seleMan2_.setSelectionSet(evaluator2_.evaluate());
91 +    }            
92 +    
93 +    count_.resize(nThetaBins_);
94 +    histogram_.resize(nThetaBins_);
95    }
96 +
97 +  pAngle::pAngle(SimInfo* info, const std::string& filename,
98 +                 const std::string& sele1, int seleOffset, int nthetabins)
99 +    : StaticAnalyser(info, filename), selectionScript1_(sele1),  
100 +      evaluator1_(info), evaluator2_(info), seleMan1_(info), seleMan2_(info),
101 +      nThetaBins_(nthetabins), seleOffset_(seleOffset),
102 +      doVect_(false), doOffset_(true), doOffset2_(false) {
103 +    
104 +    setOutputName(getPrefix(filename) + ".pAngle");
105 +    
106 +    evaluator1_.loadScriptString(sele1);
107 +    if (!evaluator1_.isDynamic()) {
108 +      seleMan1_.setSelectionSet(evaluator1_.evaluate());
109 +    }            
110 +    
111 +    count_.resize(nThetaBins_);
112 +    histogram_.resize(nThetaBins_);    
113 +  }
114 +
115 +  pAngle::pAngle(SimInfo* info, const std::string& filename,
116 +                 const std::string& sele1, int seleOffset, int seleOffset2,
117 +                 int nthetabins)
118 +    : StaticAnalyser(info, filename), selectionScript1_(sele1),  
119 +      evaluator1_(info), evaluator2_(info), seleMan1_(info), seleMan2_(info),
120 +      nThetaBins_(nthetabins), seleOffset_(seleOffset),
121 +      seleOffset2_(seleOffset2),
122 +      doVect_(false), doOffset_(true), doOffset2_(true) {
123 +    
124 +    setOutputName(getPrefix(filename) + ".pAngle");
125 +    
126 +    evaluator1_.loadScriptString(sele1);
127 +    if (!evaluator1_.isDynamic()) {
128 +      seleMan1_.setSelectionSet(evaluator1_.evaluate());
129 +    }            
130 +    
131 +    count_.resize(nThetaBins_);
132 +    histogram_.resize(nThetaBins_);    
133 +  }
134    
135    void pAngle::process() {
136      Molecule* mol;
137      RigidBody* rb;
74 –    StuntDouble* sd;
138      SimInfo::MoleculeIterator mi;
139      Molecule::RigidBodyIterator rbIter;
140 <    int i;
140 >    StuntDouble* sd1;
141 >    StuntDouble* sd2;
142 >    int ii;
143 >    int jj;
144  
145      Thermo thermo(info_);
146      DumpReader reader(info_, dumpFilename_);    
147      int nFrames = reader.getNFrames();
148 +
149      nProcessed_ = nFrames/step_;
150  
151      std::fill(histogram_.begin(), histogram_.end(), 0.0);
# Line 99 | Line 166 | namespace OpenMD {
166        
167        Vector3d CenterOfMass = thermo.getCom();      
168  
169 <      if  (evaluator_.isDynamic()) {
170 <        seleMan_.setSelectionSet(evaluator_.evaluate());
169 >      if  (evaluator1_.isDynamic()) {
170 >        seleMan1_.setSelectionSet(evaluator1_.evaluate());
171        }
172        
173 <      for (sd = seleMan_.beginSelected(i); sd != NULL;
174 <           sd = seleMan_.nextSelected(i)) {
173 >      if (doVect_) {
174 >
175          
176 <        Vector3d pos = sd->getPos();
177 <        
111 <        Vector3d r1 = CenterOfMass - pos;
112 <        // only do this if the stunt double actually has a vector associated
113 <        // with it
114 <        if (sd->isDirectional()) {
115 <          Vector3d dipole = sd->getA().getColumn(2);
116 <          RealType distance = r1.length();
176 >        for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
177 >             sd1 = seleMan1_.nextSelected(ii)) {
178            
179 <          dipole.normalize();
180 <          r1.normalize();
181 <          RealType cosangle = dot(r1, dipole);
182 <                          
183 <          int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
184 <          count_[binNo]++;
179 >          Vector3d pos = sd1->getPos();
180 >          
181 >          Vector3d r1 = CenterOfMass - pos;
182 >          // only do this if the stunt double actually has a vector associated
183 >          // with it
184 >          if (sd1->isDirectional()) {
185 >            Vector3d vec = sd1->getA().getColumn(2);
186 >            RealType distance = r1.length();
187 >            
188 >            vec.normalize();
189 >            r1.normalize();
190 >            RealType cosangle = dot(r1, vec);
191 >            
192 >            int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
193 >            count_[binNo]++;
194 >          }
195          }
196 <        
196 >      } else {
197 >        if (doOffset_) {
198 >          
199 >          for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
200 >               sd1 = seleMan1_.nextSelected(ii)) {
201 >            
202 >            // This will require careful rewriting if StaticProps is
203 >            // ever parallelized.  For an example, see
204 >            // Thermo::getTaggedAtomPairDistance
205 >            Vector3d r1;
206 >            
207 >            if (doOffset2_) {
208 >              int sd1Aind = sd1->getGlobalIndex() + seleOffset2_;
209 >              StuntDouble* sd1A = info_->getIOIndexToIntegrableObject(sd1Aind);
210 >              r1 = CenterOfMass - sd1A->getPos();              
211 >            } else {
212 >              r1 = CenterOfMass - sd1->getPos();
213 >            }
214 >
215 >            if (usePeriodicBoundaryConditions_)
216 >              currentSnapshot_->wrapVector(r1);
217 >
218 >
219 >            int sd2Index = sd1->getGlobalIndex() + seleOffset_;
220 >            sd2 = info_->getIOIndexToIntegrableObject(sd2Index);
221 >
222 >            Vector3d r2 = CenterOfMass - sd2->getPos();
223 >            if (usePeriodicBoundaryConditions_)
224 >              currentSnapshot_->wrapVector(r1);
225 >            
226 >            Vector3d rc = 0.5*(r1 + r2);
227 >            if (usePeriodicBoundaryConditions_)
228 >              currentSnapshot_->wrapVector(rc);
229 >            
230 >            Vector3d vec = r1-r2;
231 >            if (usePeriodicBoundaryConditions_)
232 >              currentSnapshot_->wrapVector(vec);
233 >            
234 >            rc.normalize();
235 >            vec.normalize();
236 >            RealType cosangle = dot(rc, vec);
237 >            int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
238 >            count_[binNo]++;
239 >          }        
240 >        } else {
241 >          
242 >          if  (evaluator2_.isDynamic()) {
243 >            seleMan2_.setSelectionSet(evaluator2_.evaluate());
244 >          }
245 >          
246 >          if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) {
247 >            sprintf( painCave.errMsg,
248 >                     "In frame %d, the number of selected StuntDoubles are\n"
249 >                     "\tnot the same in --sele1 and sele2\n", istep);
250 >            painCave.severity = OPENMD_INFO;
251 >            painCave.isFatal = 0;
252 >            simError();            
253 >          }
254 >          
255 >          for (sd1 = seleMan1_.beginSelected(ii),
256 >                 sd2 = seleMan2_.beginSelected(jj);
257 >               sd1 != NULL && sd2 != NULL;
258 >               sd1 = seleMan1_.nextSelected(ii),
259 >                 sd2 = seleMan2_.nextSelected(jj)) {
260 >            
261 >            Vector3d r1 = CenterOfMass - sd1->getPos();
262 >            if (usePeriodicBoundaryConditions_)
263 >              currentSnapshot_->wrapVector(r1);
264 >            
265 >            Vector3d r2 = CenterOfMass - sd2->getPos();
266 >            if (usePeriodicBoundaryConditions_)
267 >              currentSnapshot_->wrapVector(r1);
268 >          
269 >            Vector3d rc = 0.5*(r1 + r2);
270 >            if (usePeriodicBoundaryConditions_)
271 >              currentSnapshot_->wrapVector(rc);
272 >            
273 >            Vector3d vec = r1-r2;
274 >            if (usePeriodicBoundaryConditions_)
275 >              currentSnapshot_->wrapVector(vec);
276 >            
277 >            rc.normalize();
278 >            vec.normalize();
279 >            RealType cosangle = dot(rc, vec);
280 >            int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
281 >            count_[binNo]++;
282 >            
283 >          }
284 >        }
285        }
286      }
287 +      
288      processHistogram();
289      writeProbs();
290      
# Line 148 | Line 308 | namespace OpenMD {
308      if (rdfStream.is_open()) {
309        rdfStream << "#pAngle\n";
310        rdfStream << "#nFrames:\t" << nProcessed_ << "\n";
311 <      rdfStream << "#selection: (" << selectionScript_ << ")\n";
311 >      rdfStream << "#selection1: (" << selectionScript1_ << ")";
312 >      if (!doVect_) {
313 >        rdfStream << "\tselection2: (" << selectionScript2_ << ")";
314 >      }
315 >      rdfStream << "\n";
316        rdfStream << "#cos(theta)\tp(cos(theta))\n";
317        RealType dct = 2.0 / histogram_.size();
318        for (unsigned int i = 0; i < histogram_.size(); ++i) {
# Line 158 | Line 322 | namespace OpenMD {
322        
323      } else {
324        
325 <      sprintf(painCave.errMsg, "pAngle: unable to open %s\n", outputFilename_.c_str());
325 >      sprintf(painCave.errMsg, "pAngle: unable to open %s\n",
326 >              outputFilename_.c_str());
327        painCave.isFatal = 1;
328        simError();  
329      }

Diff Legend

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