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

Comparing trunk/src/applications/staticProps/GofAngle2.cpp (file contents):
Revision 2022 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 2023 by gezelter, Thu Oct 2 14:35:14 2014 UTC

# Line 49 | Line 49 | namespace OpenMD {
49  
50   namespace OpenMD {
51  
52 <  GofAngle2::GofAngle2(SimInfo* info, const std::string& filename, const std::string& sele1,
52 >  GofAngle2::GofAngle2(SimInfo* info, const std::string& filename,
53 >                       const std::string& sele1,
54                         const std::string& sele2, int nangleBins)
55 <    : RadialDistrFunc(info, filename, sele1, sele2), nAngleBins_(nangleBins) {
55 >    : RadialDistrFunc(info, filename, sele1, sele2), nAngleBins_(nangleBins),
56 >      evaluator3_(info),
57 >      seleMan3_(info), doSele3_(false) {
58 >    
59 >    setOutputName(getPrefix(filename) + ".gto");
60 >    
61 >    deltaCosAngle_ = 2.0 / nAngleBins_;
62 >    
63 >    histogram_.resize(nAngleBins_);
64 >    avgGofr_.resize(nAngleBins_);
65 >    for (int i = 0 ; i < nAngleBins_; ++i) {
66 >      histogram_[i].resize(nAngleBins_);
67 >      avgGofr_[i].resize(nAngleBins_);
68 >    }  
69 >  }
70  
71 <      setOutputName(getPrefix(filename) + ".gto");
71 >  GofAngle2::GofAngle2(SimInfo* info, const std::string& filename,
72 >                       const std::string& sele1,
73 >                       const std::string& sele2,
74 >                       const std::string& sele3, int nangleBins)
75 >    : RadialDistrFunc(info, filename, sele1, sele2), nAngleBins_(nangleBins),
76 >      evaluator3_(info), selectionScript3_(sele3),
77 >      seleMan3_(info), doSele3_(true) {
78 >    
79 >    setOutputName(getPrefix(filename) + ".gto");
80  
81 <      deltaCosAngle_ = 2.0 / nAngleBins_;
81 >    deltaCosAngle_ = 2.0 / nAngleBins_;
82 >    
83 >    histogram_.resize(nAngleBins_);
84 >    avgGofr_.resize(nAngleBins_);
85 >    for (int i = 0 ; i < nAngleBins_; ++i) {
86 >      histogram_[i].resize(nAngleBins_);
87 >      avgGofr_[i].resize(nAngleBins_);
88 >    }    
89 >    evaluator3_.loadScriptString(sele3);      
90 >    if (!evaluator3_.isDynamic()) {
91 >      seleMan3_.setSelectionSet(evaluator3_.evaluate());
92 >    }
93 >  }
94  
95 <      histogram_.resize(nAngleBins_);
96 <      avgGofr_.resize(nAngleBins_);
97 <      for (int i = 0 ; i < nAngleBins_; ++i) {
98 <        histogram_[i].resize(nAngleBins_);
99 <        avgGofr_[i].resize(nAngleBins_);
100 <      }    
95 >  void GofAngle2::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 GofAngle2::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 GofAngle2::preProcess() {
165  
166      for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
# Line 104 | Line 198 | namespace OpenMD {
198      MultipoleAdapter ma1 = MultipoleAdapter(atype1);
199      MultipoleAdapter ma2 = MultipoleAdapter(atype2);
200  
201 +    if (!sd1->isDirectional()) {
202 +      sprintf(painCave.errMsg,
203 +              "GofAngle2: attempted to use a non-directional object: %s\n",
204 +              sd1->getType().c_str());
205 +      painCave.isFatal = 1;
206 +      simError();  
207 +    }
208 +
209 +    if (!sd2->isDirectional()) {
210 +      sprintf(painCave.errMsg,
211 +              "GofAngle2: attempted to use a non-directional object: %s\n",
212 +              sd2->getType().c_str());
213 +      painCave.isFatal = 1;
214 +      simError();  
215 +    }
216 +
217      Vector3d dipole1, dipole2;
218      if (ma1.isDipole())        
219          dipole1 = sd1->getDipole();
# Line 131 | Line 241 | namespace OpenMD {
241      ++npairs_;
242    }
243  
244 +  void GofAngle2::collectHistogram(StuntDouble* sd1, StuntDouble* sd2,
245 +                                   StuntDouble* sd3) {
246 +
247 +    if (sd1 == sd2) {
248 +      return;
249 +    }
250 +
251 +    Vector3d p1 = sd1->getPos();
252 +    Vector3d p3 = sd3->getPos();
253 +
254 +    Vector3d c = 0.5 * (p1 + p3);
255 +    Vector3d r13 = p3 - p1;
256 +
257 +    Vector3d r12 = sd2->getPos() - c;
258 +  
259 +    if (usePeriodicBoundaryConditions_) {
260 +      currentSnapshot_->wrapVector(r12);
261 +      currentSnapshot_->wrapVector(r13);
262 +    }
263 +    r12.normalize();
264 +    r13.normalize();
265 +
266 +    if (!sd2->isDirectional()) {
267 +      sprintf(painCave.errMsg,
268 +              "GofAngle2: attempted to use a non-directional object: %s\n",
269 +              sd2->getType().c_str());
270 +      painCave.isFatal = 1;
271 +      simError();  
272 +    }
273 +
274 +    AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
275 +    MultipoleAdapter ma2 = MultipoleAdapter(atype2);
276 +
277 +    Vector3d dipole2;
278 +    if (ma2.isDipole())        
279 +        dipole2 = sd2->getDipole();
280 +    else
281 +        dipole2 = sd2->getA().transpose() * V3Z;
282 +    
283 +    dipole2.normalize();    
284 +
285 +    RealType cosAngle1 = dot(r12, r13);
286 +    RealType cosAngle2 = dot(r13, dipole2);
287 +
288 +    RealType halfBin = (nAngleBins_ - 1) * 0.5;
289 +    int angleBin1 = int(halfBin * (cosAngle1 + 1.0));
290 +    int angleBin2 = int(halfBin * (cosAngle2 + 1.0));
291 +
292 +    ++histogram_[angleBin1][angleBin2];    
293 +    ++npairs_;
294 +
295 +  }
296 +
297    void GofAngle2::writeRdf() {
298      std::ofstream rdfStream(outputFilename_.c_str());
299      if (rdfStream.is_open()) {
300        rdfStream << "#radial distribution function\n";
301        rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
302 <      rdfStream << "selection2: (" << selectionScript2_ << ")\n";
303 <      rdfStream << "#nAngleBins =" << nAngleBins_ << "deltaCosAngle = " << deltaCosAngle_ << "\n";
302 >      rdfStream << "selection2: (" << selectionScript2_ << ")";
303 >      if (doSele3_) {
304 >        rdfStream << "\tselection3: (" << selectionScript3_ << ")\n";
305 >      } else {
306 >        rdfStream << "\n";
307 >      }
308 >      rdfStream << "#nAngleBins =" << nAngleBins_ << "deltaCosAngle = "
309 >                << deltaCosAngle_ << "\n";
310        for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
311          // RealType cosAngle1 = -1.0 + (i + 0.5)*deltaCosAngle_;
312          
# Line 145 | Line 314 | namespace OpenMD {
314            // RealType cosAngle2 = -1.0 + (j + 0.5)*deltaCosAngle_;
315            rdfStream <<avgGofr_[i][j]/nProcessed_ << "\t";
316          }
148 –
317          rdfStream << "\n";
318        }
319          
320      } else {
321  
322 <      sprintf(painCave.errMsg, "GofAngle2: unable to open %s\n", outputFilename_.c_str());
322 >      sprintf(painCave.errMsg, "GofAngle2: unable to open %s\n",
323 >              outputFilename_.c_str());
324        painCave.isFatal = 1;
325        simError();  
326      }

Diff Legend

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