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 1879 by gezelter, Sun Jun 16 15:15:42 2013 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 122 | Line 220 | namespace OpenMD {
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];
252 +        
253 +      ++npairs_;
254 +    }
255 +  }
256 +
257    void GofRAngle::writeRdf() {
258      std::ofstream rdfStream(outputFilename_.c_str());
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 160 | Line 300 | namespace OpenMD {
300  
301      r12.normalize();
302  
303 <    AtomType* atype1 = static_cast<Atom*>(sd1)->getAtomType();
304 <    MultipoleAdapter ma1 = MultipoleAdapter(atype1);
305 <    Vector3d vec;
306 <    if (ma1.isDipole() )
307 <      vec = sd1->getDipole();
308 <    else
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);
316 >      
317 >      if (ma1.isDipole() )
318 >        vec = sd1->getDipole();
319 >      else
320 >        vec = sd1->getA().transpose() * V3Z;
321 >    } else {
322        vec = sd1->getA().transpose() * V3Z;
323 <    vec.normalize();    
323 >    }
324  
325 +    vec.normalize();    
326 +      
327      return dot(r12, vec);
328    }
329  
330 <  RealType GofROmega::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2) {
330 >  RealType GofRTheta::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2,
331 >                                    StuntDouble* sd3) {
332 >    Vector3d p1 = sd1->getPos();
333 >    Vector3d p3 = sd3->getPos();
334  
335 <    AtomType* atype1 = static_cast<Atom*>(sd1)->getAtomType();
336 <    AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
335 >    Vector3d c = 0.5 * (p1 + p3);
336 >    Vector3d r13 = p3 - p1;
337  
338 <    MultipoleAdapter ma1 = MultipoleAdapter(atype1);
339 <    MultipoleAdapter ma2 = MultipoleAdapter(atype2);
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 (ma1.isDipole() )
355 <      v1 = sd1->getDipole();
356 <    else
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);
365 >      if (ma1.isDipole() )
366 >        v1 = sd1->getDipole();
367 >      else
368 >        v1 = sd1->getA().transpose() * V3Z;
369 >    } else {
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 (ma2.isDipole() )
382 <      v2 = sd2->getDipole();
383 <    else
381 >    if (sd2->isAtom()) {
382 >      AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
383 >      MultipoleAdapter ma2 = MultipoleAdapter(atype2);
384 >          
385 >      if (ma2.isDipole() )
386 >        v2 = sd2->getDipole();
387 >      else
388 >        v2 = sd2->getA().transpose() * V3Z;
389 >    } else {
390        v2 = sd2->getA().transpose() * V3Z;
391 <
391 >    }
392 >      
393      v1.normalize();
394      v2.normalize();
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