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 2071 by gezelter, Sat Mar 7 21:41:51 2015 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), doVect_(true), doOffset_(false),
58 >      selectionScript1_(sele1), seleMan1_(info), seleMan2_(info),
59 >      evaluator1_(info),  evaluator2_(info),
60 >      nThetaBins_(nthetabins) {
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), doVect_(false), doOffset_(false),
77 >      selectionScript1_(sele1), selectionScript2_(sele2),
78 >      seleMan1_(info), seleMan2_(info), evaluator1_(info), evaluator2_(info),
79 >      nThetaBins_(nthetabins) {
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), doVect_(false), doOffset_(true),
100 +      doOffset2_(false), selectionScript1_(sele1),  
101 +      seleMan1_(info), seleMan2_(info), evaluator1_(info), evaluator2_(info),
102 +      seleOffset_(seleOffset),  nThetaBins_(nthetabins) {
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), doVect_(false), doOffset_(true),
119 +      doOffset2_(true), selectionScript1_(sele1),  
120 +      seleMan1_(info), seleMan2_(info), evaluator1_(info), evaluator2_(info),
121 +      seleOffset_(seleOffset), seleOffset2_(seleOffset2),
122 +      nThetaBins_(nthetabins) {
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 >            vec.normalize();
187 >            r1.normalize();
188 >            RealType cosangle = dot(r1, vec);
189 >            
190 >            int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
191 >            count_[binNo]++;
192 >          }
193          }
194 <        
194 >      } else {
195 >        if (doOffset_) {
196 >          
197 >          for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
198 >               sd1 = seleMan1_.nextSelected(ii)) {
199 >            
200 >            // This will require careful rewriting if StaticProps is
201 >            // ever parallelized.  For an example, see
202 >            // Thermo::getTaggedAtomPairDistance
203 >            Vector3d r1;
204 >            
205 >            if (doOffset2_) {
206 >              int sd1Aind = sd1->getGlobalIndex() + seleOffset2_;
207 >              StuntDouble* sd1A = info_->getIOIndexToIntegrableObject(sd1Aind);
208 >              r1 = CenterOfMass - sd1A->getPos();              
209 >            } else {
210 >              r1 = CenterOfMass - sd1->getPos();
211 >            }
212 >
213 >            if (usePeriodicBoundaryConditions_)
214 >              currentSnapshot_->wrapVector(r1);
215 >
216 >
217 >            int sd2Index = sd1->getGlobalIndex() + seleOffset_;
218 >            sd2 = info_->getIOIndexToIntegrableObject(sd2Index);
219 >
220 >            Vector3d r2 = CenterOfMass - sd2->getPos();
221 >            if (usePeriodicBoundaryConditions_)
222 >              currentSnapshot_->wrapVector(r1);
223 >            
224 >            Vector3d rc = 0.5*(r1 + r2);
225 >            if (usePeriodicBoundaryConditions_)
226 >              currentSnapshot_->wrapVector(rc);
227 >            
228 >            Vector3d vec = r1-r2;
229 >            if (usePeriodicBoundaryConditions_)
230 >              currentSnapshot_->wrapVector(vec);
231 >            
232 >            rc.normalize();
233 >            vec.normalize();
234 >            RealType cosangle = dot(rc, vec);
235 >            int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
236 >            count_[binNo]++;
237 >          }        
238 >        } else {
239 >          
240 >          if  (evaluator2_.isDynamic()) {
241 >            seleMan2_.setSelectionSet(evaluator2_.evaluate());
242 >          }
243 >          
244 >          if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) {
245 >            sprintf( painCave.errMsg,
246 >                     "In frame %d, the number of selected StuntDoubles are\n"
247 >                     "\tnot the same in --sele1 and sele2\n", istep);
248 >            painCave.severity = OPENMD_INFO;
249 >            painCave.isFatal = 0;
250 >            simError();            
251 >          }
252 >          
253 >          for (sd1 = seleMan1_.beginSelected(ii),
254 >                 sd2 = seleMan2_.beginSelected(jj);
255 >               sd1 != NULL && sd2 != NULL;
256 >               sd1 = seleMan1_.nextSelected(ii),
257 >                 sd2 = seleMan2_.nextSelected(jj)) {
258 >            
259 >            Vector3d r1 = CenterOfMass - sd1->getPos();
260 >            if (usePeriodicBoundaryConditions_)
261 >              currentSnapshot_->wrapVector(r1);
262 >            
263 >            Vector3d r2 = CenterOfMass - sd2->getPos();
264 >            if (usePeriodicBoundaryConditions_)
265 >              currentSnapshot_->wrapVector(r1);
266 >          
267 >            Vector3d rc = 0.5*(r1 + r2);
268 >            if (usePeriodicBoundaryConditions_)
269 >              currentSnapshot_->wrapVector(rc);
270 >            
271 >            Vector3d vec = r1-r2;
272 >            if (usePeriodicBoundaryConditions_)
273 >              currentSnapshot_->wrapVector(vec);
274 >            
275 >            rc.normalize();
276 >            vec.normalize();
277 >            RealType cosangle = dot(rc, vec);
278 >            int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
279 >            count_[binNo]++;
280 >            
281 >          }
282 >        }
283        }
284      }
285 +      
286      processHistogram();
287      writeProbs();
288      
# Line 148 | Line 306 | namespace OpenMD {
306      if (rdfStream.is_open()) {
307        rdfStream << "#pAngle\n";
308        rdfStream << "#nFrames:\t" << nProcessed_ << "\n";
309 <      rdfStream << "#selection: (" << selectionScript_ << ")\n";
309 >      rdfStream << "#selection1: (" << selectionScript1_ << ")";
310 >      if (!doVect_) {
311 >        rdfStream << "\tselection2: (" << selectionScript2_ << ")";
312 >      }
313 >      rdfStream << "\n";
314        rdfStream << "#cos(theta)\tp(cos(theta))\n";
315        RealType dct = 2.0 / histogram_.size();
316        for (unsigned int i = 0; i < histogram_.size(); ++i) {
# Line 158 | Line 320 | namespace OpenMD {
320        
321      } else {
322        
323 <      sprintf(painCave.errMsg, "pAngle: unable to open %s\n", outputFilename_.c_str());
323 >      sprintf(painCave.errMsg, "pAngle: unable to open %s\n",
324 >              outputFilename_.c_str());
325        painCave.isFatal = 1;
326        simError();  
327      }

Diff Legend

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