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

Comparing trunk/src/applications/staticProps/GofRAngle.cpp (file contents):
Revision 2022 by gezelter, Thu Feb 20 16:27:30 2014 UTC vs.
Revision 2023 by gezelter, Thu Oct 2 14:35:14 2014 UTC

# Line 48 | Line 48 | namespace OpenMD {
48   #include "utils/simError.h"
49  
50   namespace OpenMD {
51 +  
52 +  GofRAngle::GofRAngle(SimInfo* info, const std::string& filename,
53 +                       const std::string& sele1,
54 +                       const std::string& sele2,
55 +                       RealType len, int nrbins, int nangleBins)
56 +    : RadialDistrFunc(info, filename, sele1, sele2), len_(len),
57 +      nRBins_(nrbins), nAngleBins_(nangleBins), evaluator3_(info),
58 +      seleMan3_(info), doSele3_(false) {
59 +    
60 +    deltaR_ = len_ /(double) nRBins_;
61 +    deltaCosAngle_ = 2.0 / (double)nAngleBins_;    
62 +    histogram_.resize(nRBins_);
63 +    avgGofr_.resize(nRBins_);
64 +    for (int i = 0 ; i < nRBins_; ++i) {
65 +      histogram_[i].resize(nAngleBins_);
66 +      avgGofr_[i].resize(nAngleBins_);
67 +    }
68 +  }
69 +  
70 +  GofRAngle::GofRAngle(SimInfo* info, const std::string& filename,
71 +                       const std::string& sele1,
72 +                       const std::string& sele2,
73 +                       const std::string& sele3,
74 +                       RealType len, int nrbins, int nangleBins)
75 +    : RadialDistrFunc(info, filename, sele1, sele2), len_(len),
76 +      nRBins_(nrbins), nAngleBins_(nangleBins), selectionScript3_(sele3),
77 +      evaluator3_(info), seleMan3_(info), doSele3_(true) {
78  
79 <  GofRAngle::GofRAngle(SimInfo* info, const std::string& filename, const std::string& sele1,
80 <                       const std::string& sele2, RealType len, int nrbins, int nangleBins)
81 <    : RadialDistrFunc(info, filename, sele1, sele2), len_(len), nRBins_(nrbins), nAngleBins_(nangleBins){
79 >    deltaR_ = len_ /(double) nRBins_;
80 >    deltaCosAngle_ = 2.0 / (double)nAngleBins_;    
81 >    histogram_.resize(nRBins_);
82 >    avgGofr_.resize(nRBins_);
83 >    for (int i = 0 ; i < nRBins_; ++i) {
84 >      histogram_[i].resize(nAngleBins_);
85 >      avgGofr_[i].resize(nAngleBins_);
86 >    }
87  
88 <      deltaR_ = len_ /(double) nRBins_;
89 <      deltaCosAngle_ = 2.0 / (double)nAngleBins_;    
90 <      histogram_.resize(nRBins_);
91 <      avgGofr_.resize(nRBins_);
92 <      for (int i = 0 ; i < nRBins_; ++i) {
93 <        histogram_[i].resize(nAngleBins_);
94 <        avgGofr_[i].resize(nAngleBins_);
88 >    evaluator3_.loadScriptString(sele3);      
89 >    if (!evaluator3_.isDynamic()) {
90 >      seleMan3_.setSelectionSet(evaluator3_.evaluate());
91 >    }
92 >
93 >  }
94 >  
95 >  void GofRAngle::processNonOverlapping( SelectionManager& sman1,
96 >                                         SelectionManager& sman2) {
97 >    StuntDouble* sd1;
98 >    StuntDouble* sd2;
99 >    StuntDouble* sd3;
100 >    int i;    
101 >    int j;
102 >    int k;
103 >    
104 >    // This is the same as a non-overlapping pairwise loop structure:
105 >    // for (int i = 0;  i < ni ; ++i ) {
106 >    //   for (int j = 0; j < nj; ++j) {}
107 >    // }
108 >
109 >    if (doSele3_) {
110 >      if  (evaluator3_.isDynamic()) {
111 >        seleMan3_.setSelectionSet(evaluator3_.evaluate());
112        }
113 +      if (sman1.getSelectionCount() != seleMan3_.getSelectionCount() ) {
114 +        RadialDistrFunc::processNonOverlapping( sman1, sman2 );
115 +      }
116 +
117 +      for (sd1 = sman1.beginSelected(i), sd3 = seleMan3_.beginSelected(k);
118 +           sd1 != NULL && sd3 != NULL;
119 +           sd1 = sman1.nextSelected(i), sd3 = seleMan3_.nextSelected(k)) {
120 +        for (sd2 = sman2.beginSelected(j); sd2 != NULL;
121 +             sd2 = sman2.nextSelected(j)) {
122 +          collectHistogram(sd1, sd2, sd3);
123 +        }
124 +      }
125 +    } else {
126 +      RadialDistrFunc::processNonOverlapping( sman1, sman2 );
127      }
128 +  }
129  
130 +  void GofRAngle::processOverlapping( SelectionManager& sman) {
131 +    StuntDouble* sd1;
132 +    StuntDouble* sd2;
133 +    StuntDouble* sd3;
134 +    int i;    
135 +    int j;
136 +    int k;
137  
138 +    // This is the same as a pairwise loop structure:
139 +    // for (int i = 0;  i < n-1 ; ++i ) {
140 +    //   for (int j = i + 1; j < n; ++j) {}
141 +    // }
142 +    
143 +    if (doSele3_) {
144 +      if  (evaluator3_.isDynamic()) {
145 +        seleMan3_.setSelectionSet(evaluator3_.evaluate());
146 +      }
147 +      if (sman.getSelectionCount() != seleMan3_.getSelectionCount() ) {
148 +        RadialDistrFunc::processOverlapping( sman);
149 +      }
150 +      for (sd1 = sman.beginSelected(i), sd3 = seleMan3_.beginSelected(k);
151 +           sd1 != NULL && sd3 != NULL;
152 +           sd1 = sman.nextSelected(i), sd3 = seleMan3_.nextSelected(k)) {
153 +        for (j  = i, sd2 = sman.nextSelected(j); sd2 != NULL;
154 +             sd2 = sman.nextSelected(j)) {
155 +          collectHistogram(sd1, sd2, sd3);
156 +        }            
157 +      }
158 +    } else {
159 +      RadialDistrFunc::processOverlapping( sman);
160 +    }    
161 +  }
162 +  
163 +
164    void GofRAngle::preProcess() {
165      for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
166        std::fill(avgGofr_[i].begin(), avgGofr_[i].end(), 0);
# Line 87 | Line 184 | namespace OpenMD {
184  
185        RealType rLower = i * deltaR_;
186        RealType rUpper = rLower + deltaR_;
187 <      RealType volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower );
187 >      RealType volSlice = ( rUpper * rUpper * rUpper ) -
188 >        ( rLower * rLower * rLower );
189        RealType nIdeal = volSlice * pairConstant;
190  
191        for (unsigned int j = 0; j < histogram_[i].size(); ++j){
# Line 114 | Line 212 | namespace OpenMD {
212      if (distance <= len_) {
213  
214        RealType cosAngle = evaluateAngle(sd1, sd2);
215 +      RealType halfBin = (nAngleBins_ - 1) * 0.5;
216 +      int whichThetaBin = int(halfBin * (cosAngle + 1.0));
217 +      ++histogram_[whichRBin][whichThetaBin];
218 +        
219 +      ++npairs_;
220 +    }
221 +  }
222 +
223 +  void GofRAngle::collectHistogram(StuntDouble* sd1, StuntDouble* sd2,
224 +                                   StuntDouble* sd3) {
225 +
226 +    if (sd1 == sd2) {
227 +      return;
228 +    }
229 +
230 +    Vector3d p1 = sd1->getPos();
231 +    Vector3d p3 = sd3->getPos();
232 +
233 +    Vector3d c = 0.5 * (p1 + p3);
234 +    Vector3d r13 = p3 - p1;
235 +
236 +    Vector3d r12 = sd2->getPos() - c;
237 +  
238 +    if (usePeriodicBoundaryConditions_) {
239 +      currentSnapshot_->wrapVector(r12);
240 +      currentSnapshot_->wrapVector(r13);
241 +    }
242 +    
243 +    RealType distance = r12.length();
244 +    int whichRBin = int(distance / deltaR_);
245 +
246 +    if (distance <= len_) {
247 +
248 +      RealType cosAngle = evaluateAngle(sd1, sd2, sd3);
249        RealType halfBin = (nAngleBins_ - 1) * 0.5;
250        int whichThetaBin = int(halfBin * (cosAngle + 1.0));
251        ++histogram_[whichRBin][whichThetaBin];
# Line 127 | Line 259 | namespace OpenMD {
259      if (rdfStream.is_open()) {
260        rdfStream << "#radial distribution function\n";
261        rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
262 <      rdfStream << "selection2: (" << selectionScript2_ << ")\n";
263 <      rdfStream << "#nRBins = " << nRBins_ << "\t maxLen = " << len_ << "deltaR = " << deltaR_ <<"\n";
264 <      rdfStream << "#nAngleBins =" << nAngleBins_ << "deltaCosAngle = " << deltaCosAngle_ << "\n";
262 >      rdfStream <<  "selection2: (" << selectionScript2_ << ")";
263 >      if (doSele3_) {
264 >        rdfStream << "\tselection3: (" << selectionScript3_ << ")\n";
265 >      } else {
266 >        rdfStream << "\n";
267 >      }
268 >      rdfStream << "#nRBins = " << nRBins_ << "\tmaxLen = "
269 >                << len_ << "\tdeltaR = " << deltaR_ <<"\n";
270 >      rdfStream << "#nAngleBins =" << nAngleBins_ << "\tdeltaCosAngle = "
271 >                << deltaCosAngle_ << "\n";
272        for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
273          // RealType r = deltaR_ * (i + 0.5);
274  
# Line 142 | Line 281 | namespace OpenMD {
281        }
282          
283      } else {
284 <      sprintf(painCave.errMsg, "GofRAngle: unable to open %s\n", outputFilename_.c_str());
284 >      sprintf(painCave.errMsg, "GofRAngle: unable to open %s\n",
285 >              outputFilename_.c_str());
286        painCave.isFatal = 1;
287        simError();  
288      }
# Line 162 | Line 302 | namespace OpenMD {
302  
303      Vector3d vec;    
304      
305 +    if (!sd1->isDirectional()) {
306 +      sprintf(painCave.errMsg,
307 +              "GofRTheta: attempted to use a non-directional object: %s\n",
308 +              sd1->getType().c_str());
309 +      painCave.isFatal = 1;
310 +      simError();  
311 +    }
312 +
313      if (sd1->isAtom()) {
314        AtomType* atype1 = static_cast<Atom*>(sd1)->getAtomType();
315        MultipoleAdapter ma1 = MultipoleAdapter(atype1);
# Line 179 | Line 327 | namespace OpenMD {
327      return dot(r12, vec);
328    }
329  
330 +  RealType GofRTheta::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2,
331 +                                    StuntDouble* sd3) {
332 +    Vector3d p1 = sd1->getPos();
333 +    Vector3d p3 = sd3->getPos();
334 +
335 +    Vector3d c = 0.5 * (p1 + p3);
336 +    Vector3d r13 = p3 - p1;
337 +
338 +    Vector3d r12 = sd2->getPos() - c;
339 +  
340 +    if (usePeriodicBoundaryConditions_) {
341 +      currentSnapshot_->wrapVector(r12);
342 +      currentSnapshot_->wrapVector(r13);
343 +    }
344 +
345 +    r12.normalize();
346 +    r13.normalize();
347 +          
348 +    return dot(r12, r13);
349 +  }
350 +
351    RealType GofROmega::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2) {
352      Vector3d v1, v2;
353  
354 +    if (!sd1->isDirectional()) {
355 +      sprintf(painCave.errMsg,
356 +              "GofROmega: attempted to use a non-directional object: %s\n",
357 +              sd1->getType().c_str());
358 +      painCave.isFatal = 1;
359 +      simError();  
360 +    }
361 +
362      if (sd1->isAtom()){
363        AtomType* atype1 = static_cast<Atom*>(sd1)->getAtomType();
364        MultipoleAdapter ma1 = MultipoleAdapter(atype1);
# Line 193 | Line 370 | namespace OpenMD {
370        v1 = sd1->getA().transpose() * V3Z;
371      }
372      
373 +    if (!sd2->isDirectional()) {
374 +      sprintf(painCave.errMsg,
375 +              "GofROmega attempted to use a non-directional object: %s\n",
376 +              sd2->getType().c_str());
377 +      painCave.isFatal = 1;
378 +      simError();  
379 +    }
380 +
381      if (sd2->isAtom()) {
382        AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
383        MultipoleAdapter ma2 = MultipoleAdapter(atype2);
# Line 210 | Line 395 | namespace OpenMD {
395      return dot(v1, v2);
396    }
397  
398 +  RealType GofROmega::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2,
399 +                                    StuntDouble* sd3) {
400  
401 +    Vector3d v1;
402 +    Vector3d v2;
403 +    
404 +    v1 = sd3->getPos() - sd1->getPos();
405 +    if (usePeriodicBoundaryConditions_)
406 +      currentSnapshot_->wrapVector(v1);
407 +    
408 +    if (!sd2->isDirectional()) {
409 +      sprintf(painCave.errMsg,
410 +              "GofROmega: attempted to use a non-directional object: %s\n",
411 +              sd2->getType().c_str());
412 +      painCave.isFatal = 1;
413 +      simError();  
414 +    }
415 +
416 +    if (sd2->isAtom()) {
417 +      AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
418 +      MultipoleAdapter ma2 = MultipoleAdapter(atype2);
419 +          
420 +      if (ma2.isDipole() )
421 +        v2 = sd2->getDipole();
422 +      else
423 +        v2 = sd2->getA().transpose() * V3Z;
424 +    } else {
425 +      v2 = sd2->getA().transpose() * V3Z;
426 +    }
427 +      
428 +    v1.normalize();
429 +    v2.normalize();
430 +    return dot(v1, v2);
431 +  }
432   }
433  
434  

Diff Legend

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