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 1796 by gezelter, Mon Sep 10 18:38:44 2012 UTC vs.
Revision 2023 by gezelter, Thu Oct 2 14:35:14 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   * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
# Line 43 | Line 43
43   #include <algorithm>
44   #include <fstream>
45   #include "applications/staticProps/GofAngle2.hpp"
46 + #include "primitives/Atom.hpp"
47 + #include "types/MultipoleAdapter.hpp"
48   #include "utils/simError.h"
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 96 | Line 192 | namespace OpenMD {
192      Vector3d r12 = pos1 - pos2;
193      if (usePeriodicBoundaryConditions_)
194        currentSnapshot_->wrapVector(r12);
195 <    Vector3d dipole1 = sd1->getElectroFrame().getColumn(2);
196 <    Vector3d dipole2 = sd2->getElectroFrame().getColumn(2);
195 >
196 >    AtomType* atype1 = static_cast<Atom*>(sd1)->getAtomType();
197 >    AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
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();
220 >    else
221 >        dipole1 = sd1->getA().transpose() * V3Z;
222 >
223 >    if (ma2.isDipole())        
224 >        dipole2 = sd2->getDipole();
225 >    else
226 >        dipole2 = sd2->getA().transpose() * V3Z;
227      
228      r12.normalize();
229      dipole1.normalize();    
# Line 115 | 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 129 | Line 314 | namespace OpenMD {
314            // RealType cosAngle2 = -1.0 + (j + 0.5)*deltaCosAngle_;
315            rdfStream <<avgGofr_[i][j]/nProcessed_ << "\t";
316          }
132 –
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