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 1979 by gezelter, Sat Apr 5 20:56:01 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),      
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;
74 –    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 99 | Line 146 | namespace OpenMD {
146        
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 <      for (sd = seleMan_.beginSelected(i); sd != NULL;
107 <           sd = seleMan_.nextSelected(i)) {
153 >      if (doVect_) {
154          
155 <        Vector3d pos = sd->getPos();
156 <        
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();
155 >        for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
156 >             sd1 = seleMan1_.nextSelected(ii)) {
157            
158 <          dipole.normalize();
159 <          r1.normalize();
160 <          RealType cosangle = dot(r1, dipole);
161 <                          
162 <          int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
163 <          count_[binNo]++;
158 >          Vector3d pos = sd1->getPos();
159 >          
160 >          Vector3d r1 = CenterOfMass - pos;
161 >          // only do this if the stunt double actually has a vector associated
162 >          // with it
163 >          if (sd1->isDirectional()) {
164 >            Vector3d vec = sd1->getA().getColumn(2);
165 >            RealType distance = r1.length();
166 >            
167 >            vec.normalize();
168 >            r1.normalize();
169 >            RealType cosangle = dot(r1, vec);
170 >            
171 >            int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
172 >            count_[binNo]++;
173 >          }
174          }
175 <        
175 >      } else {
176 >        if (doOffset_) {
177 >          
178 >          for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
179 >               sd1 = seleMan1_.nextSelected(ii)) {
180 >            
181 >            // This will require careful rewriting if StaticProps is
182 >            // ever parallelized.  For an example, see
183 >            // Thermo::getTaggedAtomPairDistance
184 >            
185 >            int sd2Index = sd1->getGlobalIndex() + seleOffset_;
186 >            sd2 = info_->getIOIndexToIntegrableObject(sd2Index);
187 >            
188 >            Vector3d r1 = CenterOfMass - sd1->getPos();
189 >            if (usePeriodicBoundaryConditions_)
190 >              currentSnapshot_->wrapVector(r1);
191 >            
192 >            Vector3d r2 = CenterOfMass - sd2->getPos();
193 >            if (usePeriodicBoundaryConditions_)
194 >              currentSnapshot_->wrapVector(r1);
195 >            
196 >            Vector3d rc = 0.5*(r1 + r2);
197 >            if (usePeriodicBoundaryConditions_)
198 >              currentSnapshot_->wrapVector(rc);
199 >            
200 >            Vector3d vec = r1-r2;
201 >            if (usePeriodicBoundaryConditions_)
202 >              currentSnapshot_->wrapVector(vec);
203 >            
204 >            rc.normalize();
205 >            vec.normalize();
206 >            RealType cosangle = dot(rc, vec);
207 >            int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
208 >            count_[binNo]++;
209 >          }        
210 >        } else {
211 >          
212 >          if  (evaluator2_.isDynamic()) {
213 >            seleMan2_.setSelectionSet(evaluator2_.evaluate());
214 >          }
215 >          
216 >          if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) {
217 >            sprintf( painCave.errMsg,
218 >                     "In frame %d, the number of selected StuntDoubles are\n"
219 >                     "\tnot the same in --sele1 and sele2\n", istep);
220 >            painCave.severity = OPENMD_INFO;
221 >            painCave.isFatal = 0;
222 >            simError();            
223 >          }
224 >          
225 >          for (sd1 = seleMan1_.beginSelected(ii),
226 >                 sd2 = seleMan2_.beginSelected(jj);
227 >               sd1 != NULL && sd2 != NULL;
228 >               sd1 = seleMan1_.nextSelected(ii),
229 >                 sd2 = seleMan2_.nextSelected(jj)) {
230 >            
231 >            Vector3d r1 = CenterOfMass - sd1->getPos();
232 >            if (usePeriodicBoundaryConditions_)
233 >              currentSnapshot_->wrapVector(r1);
234 >            
235 >            Vector3d r2 = CenterOfMass - sd2->getPos();
236 >            if (usePeriodicBoundaryConditions_)
237 >              currentSnapshot_->wrapVector(r1);
238 >          
239 >            Vector3d rc = 0.5*(r1 + r2);
240 >            if (usePeriodicBoundaryConditions_)
241 >              currentSnapshot_->wrapVector(rc);
242 >            
243 >            Vector3d vec = r1-r2;
244 >            if (usePeriodicBoundaryConditions_)
245 >              currentSnapshot_->wrapVector(vec);
246 >            
247 >            rc.normalize();
248 >            vec.normalize();
249 >            RealType cosangle = dot(rc, vec);
250 >            int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
251 >            count_[binNo]++;
252 >            
253 >          }
254 >        }
255        }
256      }
257 +      
258      processHistogram();
259      writeProbs();
260      
# Line 148 | Line 278 | namespace OpenMD {
278      if (rdfStream.is_open()) {
279        rdfStream << "#pAngle\n";
280        rdfStream << "#nFrames:\t" << nProcessed_ << "\n";
281 <      rdfStream << "#selection: (" << selectionScript_ << ")\n";
281 >      rdfStream << "#selection1: (" << selectionScript1_ << ")";
282 >      if (!doVect_) {
283 >        rdfStream << "\tselection2: (" << selectionScript2_ << ")";
284 >      }
285 >      rdfStream << "\n";
286        rdfStream << "#cos(theta)\tp(cos(theta))\n";
287        RealType dct = 2.0 / histogram_.size();
288        for (unsigned int i = 0; i < histogram_.size(); ++i) {
# Line 158 | Line 292 | namespace OpenMD {
292        
293      } else {
294        
295 <      sprintf(painCave.errMsg, "pAngle: unable to open %s\n", outputFilename_.c_str());
295 >      sprintf(painCave.errMsg, "pAngle: unable to open %s\n",
296 >              outputFilename_.c_str());
297        painCave.isFatal = 1;
298        simError();  
299      }

Diff Legend

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