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 1522 by kstocke1, Fri Nov 19 20:26:36 2010 UTC vs.
Revision 1988 by gezelter, Fri Apr 18 19:49:54 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) {
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    void pAngle::process() {
116      Molecule* mol;
117      RigidBody* rb;
72 –    StuntDouble* sd;
118      SimInfo::MoleculeIterator mi;
119      Molecule::RigidBodyIterator rbIter;
120 <    int i;
120 >    StuntDouble* sd1;
121 >    StuntDouble* sd2;
122 >    int ii;
123 >    int jj;
124  
125 +    Thermo thermo(info_);
126      DumpReader reader(info_, dumpFilename_);    
127      int nFrames = reader.getNFrames();
128 +
129      nProcessed_ = nFrames/step_;
130  
131      std::fill(histogram_.begin(), histogram_.end(), 0.0);
# Line 94 | Line 144 | namespace OpenMD {
144          }
145        }
146        
147 <      Vector3d CenterOfMass = info_->getCom();      
147 >      Vector3d CenterOfMass = thermo.getCom();      
148  
149 <      if  (evaluator_.isDynamic()) {
150 <        seleMan_.setSelectionSet(evaluator_.evaluate());
149 >      if  (evaluator1_.isDynamic()) {
150 >        seleMan1_.setSelectionSet(evaluator1_.evaluate());
151        }
152        
153 +      if (doVect_) {
154  
104 –      int runningTot = 0;
105 –      for (sd = seleMan_.beginSelected(i); sd != NULL;
106 –           sd = seleMan_.nextSelected(i)) {
107 –      
108 –        Vector3d pos = sd->getPos();
155          
156 <        Vector3d r1 = CenterOfMass - pos;
157 <        // only do this if the stunt double actually has a vector associated
112 <        // with it
113 <        if (sd->isDirectional()) {
114 <          Vector3d dipole = sd->getA().getColumn(2);
115 <          RealType distance = r1.length();
156 >        for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
157 >             sd1 = seleMan1_.nextSelected(ii)) {
158            
159 <          dipole.normalize();
160 <          r1.normalize();
161 <          RealType cosangle = dot(r1, dipole);
162 <                          
163 <          int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
164 <          count_[binNo]++;
159 >          Vector3d pos = sd1->getPos();
160 >          
161 >          Vector3d r1 = CenterOfMass - pos;
162 >          // only do this if the stunt double actually has a vector associated
163 >          // with it
164 >          if (sd1->isDirectional()) {
165 >            Vector3d vec = sd1->getA().getColumn(2);
166 >            RealType distance = r1.length();
167 >            
168 >            vec.normalize();
169 >            r1.normalize();
170 >            RealType cosangle = dot(r1, vec);
171 >            
172 >            int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
173 >            count_[binNo]++;
174 >          }
175          }
176 <        
176 >      } else {
177 >        if (doOffset_) {
178 >          
179 >          for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
180 >               sd1 = seleMan1_.nextSelected(ii)) {
181 >            
182 >            // This will require careful rewriting if StaticProps is
183 >            // ever parallelized.  For an example, see
184 >            // Thermo::getTaggedAtomPairDistance
185 >            
186 >            int sd2Index = sd1->getGlobalIndex() + seleOffset_;
187 >            sd2 = info_->getIOIndexToIntegrableObject(sd2Index);
188 >            
189 >            Vector3d r1 = CenterOfMass - sd1->getPos();
190 >            if (usePeriodicBoundaryConditions_)
191 >              currentSnapshot_->wrapVector(r1);
192 >            
193 >            Vector3d r2 = CenterOfMass - sd2->getPos();
194 >            if (usePeriodicBoundaryConditions_)
195 >              currentSnapshot_->wrapVector(r1);
196 >            
197 >            Vector3d rc = 0.5*(r1 + r2);
198 >            if (usePeriodicBoundaryConditions_)
199 >              currentSnapshot_->wrapVector(rc);
200 >            
201 >            Vector3d vec = r1-r2;
202 >            if (usePeriodicBoundaryConditions_)
203 >              currentSnapshot_->wrapVector(vec);
204 >            
205 >            rc.normalize();
206 >            vec.normalize();
207 >            RealType cosangle = dot(rc, vec);
208 >            int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
209 >            count_[binNo]++;
210 >          }        
211 >        } else {
212 >          
213 >          if  (evaluator2_.isDynamic()) {
214 >            seleMan2_.setSelectionSet(evaluator2_.evaluate());
215 >          }
216 >          
217 >          if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) {
218 >            sprintf( painCave.errMsg,
219 >                     "In frame %d, the number of selected StuntDoubles are\n"
220 >                     "\tnot the same in --sele1 and sele2\n", istep);
221 >            painCave.severity = OPENMD_INFO;
222 >            painCave.isFatal = 0;
223 >            simError();            
224 >          }
225 >          
226 >          for (sd1 = seleMan1_.beginSelected(ii),
227 >                 sd2 = seleMan2_.beginSelected(jj);
228 >               sd1 != NULL && sd2 != NULL;
229 >               sd1 = seleMan1_.nextSelected(ii),
230 >                 sd2 = seleMan2_.nextSelected(jj)) {
231 >            
232 >            Vector3d r1 = CenterOfMass - sd1->getPos();
233 >            if (usePeriodicBoundaryConditions_)
234 >              currentSnapshot_->wrapVector(r1);
235 >            
236 >            Vector3d r2 = CenterOfMass - sd2->getPos();
237 >            if (usePeriodicBoundaryConditions_)
238 >              currentSnapshot_->wrapVector(r1);
239 >          
240 >            Vector3d rc = 0.5*(r1 + r2);
241 >            if (usePeriodicBoundaryConditions_)
242 >              currentSnapshot_->wrapVector(rc);
243 >            
244 >            Vector3d vec = r1-r2;
245 >            if (usePeriodicBoundaryConditions_)
246 >              currentSnapshot_->wrapVector(vec);
247 >            
248 >            rc.normalize();
249 >            vec.normalize();
250 >            RealType cosangle = dot(rc, vec);
251 >            int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
252 >            count_[binNo]++;
253 >            
254 >          }
255 >        }
256        }
257      }
258 +      
259      processHistogram();
260      writeProbs();
261      
# Line 132 | Line 264 | namespace OpenMD {
264    void pAngle::processHistogram() {
265      
266      int atot = 0;
267 <    for(int i = 0; i < count_.size(); ++i)
267 >    for(unsigned int i = 0; i < count_.size(); ++i)
268        atot += count_[i];
269      
270 <    for(int i = 0; i < count_.size(); ++i) {
271 <      histogram_[i] = double(count_[i]) / double(atot);
270 >    for(unsigned int i = 0; i < count_.size(); ++i) {
271 >      histogram_[i] = double(count_[i] / double(atot));
272      }    
273    }
274    
# Line 147 | Line 279 | namespace OpenMD {
279      if (rdfStream.is_open()) {
280        rdfStream << "#pAngle\n";
281        rdfStream << "#nFrames:\t" << nProcessed_ << "\n";
282 <      rdfStream << "#selection: (" << selectionScript_ << ")\n";
282 >      rdfStream << "#selection1: (" << selectionScript1_ << ")";
283 >      if (!doVect_) {
284 >        rdfStream << "\tselection2: (" << selectionScript2_ << ")";
285 >      }
286 >      rdfStream << "\n";
287        rdfStream << "#cos(theta)\tp(cos(theta))\n";
288        RealType dct = 2.0 / histogram_.size();
289 <      for (int i = 0; i < histogram_.size(); ++i) {
289 >      for (unsigned int i = 0; i < histogram_.size(); ++i) {
290          RealType ct = -1.0 + (2.0 * i + 1) / (histogram_.size());
291          rdfStream << ct << "\t" << histogram_[i]/dct << "\n";
292        }
293        
294      } else {
295        
296 <      sprintf(painCave.errMsg, "pAngle: unable to open %s\n", outputFilename_.c_str());
296 >      sprintf(painCave.errMsg, "pAngle: unable to open %s\n",
297 >              outputFilename_.c_str());
298        painCave.isFatal = 1;
299        simError();  
300      }

Diff Legend

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