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 1442 by gezelter, Mon May 10 17:28:26 2010 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).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
40 < *
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   */
42  
43   /* Calculates Rho(theta) */
44  
45   #include <algorithm>
46 <  #include <fstream>
46 > #include <fstream>
47   #include "applications/staticProps/pAngle.hpp"
48   #include "utils/simError.h"
49   #include "io/DumpReader.hpp"
50   #include "primitives/Molecule.hpp"
51 + #include "brains/Thermo.hpp"
52 +
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;
72 –    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 94 | Line 164 | namespace OpenMD {
164          }
165        }
166        
167 <      Vector3d CenterOfMass = info_->getCom();      
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 <        
108 <        Vector3d r1 = CenterOfMass - pos;
109 <        // only do this if the stunt double actually has a vector associated
110 <        // with it
111 <        if (sd->isDirectional()) {
112 <          Vector3d dipole = sd->getA().getColumn(2);
113 <          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 130 | Line 293 | namespace OpenMD {
293    void pAngle::processHistogram() {
294      
295      int atot = 0;
296 <    for(int i = 0; i < count_.size(); ++i)
296 >    for(unsigned int i = 0; i < count_.size(); ++i)
297        atot += count_[i];
298      
299 <    for(int i = 0; i < count_.size(); ++i) {
300 <      histogram_[i] = double(count_[i] / atot);
299 >    for(unsigned int i = 0; i < count_.size(); ++i) {
300 >      histogram_[i] = double(count_[i] / double(atot));
301      }    
302    }
303    
# Line 145 | 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 <      for (int i = 0; i < histogram_.size(); ++i) {
318 <        RealType ct = -1.0 + i / histogram_.size();
319 <        rdfStream << ct << "\t" << histogram_[i] << "\n";
317 >      RealType dct = 2.0 / histogram_.size();
318 >      for (unsigned int i = 0; i < histogram_.size(); ++i) {
319 >        RealType ct = -1.0 + (2.0 * i + 1) / (histogram_.size());
320 >        rdfStream << ct << "\t" << histogram_[i]/dct << "\n";
321        }
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