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 2071 by gezelter, Sat Mar 7 21:41:51 2015 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,
51 >  
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 >    doSele3_(false), seleMan3_(info), evaluator3_(info) {
57 >    
58 >    setOutputName(getPrefix(filename) + ".gto");
59 >    
60 >    deltaCosAngle_ = 2.0 / nAngleBins_;
61 >    
62 >    histogram_.resize(nAngleBins_);
63 >    avgGofr_.resize(nAngleBins_);
64 >    for (int i = 0 ; i < nAngleBins_; ++i) {
65 >      histogram_[i].resize(nAngleBins_);
66 >      avgGofr_[i].resize(nAngleBins_);
67 >    }  
68 >  }
69 >  
70 >  GofAngle2::GofAngle2(SimInfo* info, const std::string& filename,
71 >                       const std::string& sele1,
72 >                       const std::string& sele2,
73 >                       const std::string& sele3, int nangleBins)
74 >  : RadialDistrFunc(info, filename, sele1, sele2), nAngleBins_(nangleBins),
75 >    doSele3_(true), seleMan3_(info), evaluator3_(info),
76 >    selectionScript3_(sele3) {
77 >    
78 >    setOutputName(getPrefix(filename) + ".gto");
79 >    
80 >    deltaCosAngle_ = 2.0 / nAngleBins_;
81 >    
82 >    histogram_.resize(nAngleBins_);
83 >    avgGofr_.resize(nAngleBins_);
84 >    for (int i = 0 ; i < nAngleBins_; ++i) {
85 >      histogram_[i].resize(nAngleBins_);
86 >      avgGofr_[i].resize(nAngleBins_);
87 >    }    
88 >    evaluator3_.loadScriptString(sele3);      
89 >    if (!evaluator3_.isDynamic()) {
90 >      seleMan3_.setSelectionSet(evaluator3_.evaluate());
91 >    }
92 >  }
93  
94 <      setOutputName(getPrefix(filename) + ".gto");
94 >  void GofAngle2::processNonOverlapping( SelectionManager& sman1,
95 >                                         SelectionManager& sman2) {
96 >    StuntDouble* sd1;
97 >    StuntDouble* sd2;
98 >    StuntDouble* sd3;
99 >    int i;    
100 >    int j;
101 >    int k;
102 >    
103 >    // This is the same as a non-overlapping pairwise loop structure:
104 >    // for (int i = 0;  i < ni ; ++i ) {
105 >    //   for (int j = 0; j < nj; ++j) {}
106 >    // }
107  
108 <      deltaCosAngle_ = 2.0 / nAngleBins_;
108 >    if (doSele3_) {
109 >      if  (evaluator3_.isDynamic()) {
110 >        seleMan3_.setSelectionSet(evaluator3_.evaluate());
111 >      }
112 >      if (sman1.getSelectionCount() != seleMan3_.getSelectionCount() ) {
113 >        RadialDistrFunc::processNonOverlapping( sman1, sman2 );
114 >      }
115  
116 <      histogram_.resize(nAngleBins_);
117 <      avgGofr_.resize(nAngleBins_);
118 <      for (int i = 0 ; i < nAngleBins_; ++i) {
119 <        histogram_[i].resize(nAngleBins_);
120 <        avgGofr_[i].resize(nAngleBins_);
121 <      }    
122 <
116 >      for (sd1 = sman1.beginSelected(i), sd3 = seleMan3_.beginSelected(k);
117 >           sd1 != NULL && sd3 != NULL;
118 >           sd1 = sman1.nextSelected(i), sd3 = seleMan3_.nextSelected(k)) {
119 >        for (sd2 = sman2.beginSelected(j); sd2 != NULL;
120 >             sd2 = sman2.nextSelected(j)) {
121 >          collectHistogram(sd1, sd2, sd3);
122 >        }
123 >      }
124 >    } else {
125 >      RadialDistrFunc::processNonOverlapping( sman1, sman2 );
126      }
127 +  }
128  
129 +  void GofAngle2::processOverlapping( SelectionManager& sman) {
130 +    StuntDouble* sd1;
131 +    StuntDouble* sd2;
132 +    StuntDouble* sd3;
133 +    int i;    
134 +    int j;
135 +    int k;
136  
137 +    // This is the same as a pairwise loop structure:
138 +    // for (int i = 0;  i < n-1 ; ++i ) {
139 +    //   for (int j = i + 1; j < n; ++j) {}
140 +    // }
141 +    
142 +    if (doSele3_) {
143 +      if  (evaluator3_.isDynamic()) {
144 +        seleMan3_.setSelectionSet(evaluator3_.evaluate());
145 +      }
146 +      if (sman.getSelectionCount() != seleMan3_.getSelectionCount() ) {
147 +        RadialDistrFunc::processOverlapping( sman);
148 +      }
149 +      for (sd1 = sman.beginSelected(i), sd3 = seleMan3_.beginSelected(k);
150 +           sd1 != NULL && sd3 != NULL;
151 +           sd1 = sman.nextSelected(i), sd3 = seleMan3_.nextSelected(k)) {
152 +        for (j  = i, sd2 = sman.nextSelected(j); sd2 != NULL;
153 +             sd2 = sman.nextSelected(j)) {
154 +          collectHistogram(sd1, sd2, sd3);
155 +        }            
156 +      }
157 +    } else {
158 +      RadialDistrFunc::processOverlapping( sman);
159 +    }    
160 +  }
161 +
162 +
163    void GofAngle2::preProcess() {
164  
165      for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
# Line 96 | Line 191 | namespace OpenMD {
191      Vector3d r12 = pos1 - pos2;
192      if (usePeriodicBoundaryConditions_)
193        currentSnapshot_->wrapVector(r12);
194 <    Vector3d dipole1 = sd1->getElectroFrame().getColumn(2);
195 <    Vector3d dipole2 = sd2->getElectroFrame().getColumn(2);
194 >
195 >    AtomType* atype1 = static_cast<Atom*>(sd1)->getAtomType();
196 >    AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
197 >    MultipoleAdapter ma1 = MultipoleAdapter(atype1);
198 >    MultipoleAdapter ma2 = MultipoleAdapter(atype2);
199 >
200 >    if (!sd1->isDirectional()) {
201 >      sprintf(painCave.errMsg,
202 >              "GofAngle2: attempted to use a non-directional object: %s\n",
203 >              sd1->getType().c_str());
204 >      painCave.isFatal = 1;
205 >      simError();  
206 >    }
207 >
208 >    if (!sd2->isDirectional()) {
209 >      sprintf(painCave.errMsg,
210 >              "GofAngle2: attempted to use a non-directional object: %s\n",
211 >              sd2->getType().c_str());
212 >      painCave.isFatal = 1;
213 >      simError();  
214 >    }
215 >
216 >    Vector3d dipole1, dipole2;
217 >    if (ma1.isDipole())        
218 >        dipole1 = sd1->getDipole();
219 >    else
220 >        dipole1 = sd1->getA().transpose() * V3Z;
221 >
222 >    if (ma2.isDipole())        
223 >        dipole2 = sd2->getDipole();
224 >    else
225 >        dipole2 = sd2->getA().transpose() * V3Z;
226      
227      r12.normalize();
228      dipole1.normalize();    
# Line 115 | Line 240 | namespace OpenMD {
240      ++npairs_;
241    }
242  
243 +  void GofAngle2::collectHistogram(StuntDouble* sd1, StuntDouble* sd2,
244 +                                   StuntDouble* sd3) {
245 +
246 +    if (sd1 == sd2) {
247 +      return;
248 +    }
249 +
250 +    Vector3d p1 = sd1->getPos();
251 +    Vector3d p3 = sd3->getPos();
252 +
253 +    Vector3d c = 0.5 * (p1 + p3);
254 +    Vector3d r13 = p3 - p1;
255 +
256 +    Vector3d r12 = sd2->getPos() - c;
257 +  
258 +    if (usePeriodicBoundaryConditions_) {
259 +      currentSnapshot_->wrapVector(r12);
260 +      currentSnapshot_->wrapVector(r13);
261 +    }
262 +    r12.normalize();
263 +    r13.normalize();
264 +
265 +    if (!sd2->isDirectional()) {
266 +      sprintf(painCave.errMsg,
267 +              "GofAngle2: attempted to use a non-directional object: %s\n",
268 +              sd2->getType().c_str());
269 +      painCave.isFatal = 1;
270 +      simError();  
271 +    }
272 +
273 +    AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
274 +    MultipoleAdapter ma2 = MultipoleAdapter(atype2);
275 +
276 +    Vector3d dipole2;
277 +    if (ma2.isDipole())        
278 +        dipole2 = sd2->getDipole();
279 +    else
280 +        dipole2 = sd2->getA().transpose() * V3Z;
281 +    
282 +    dipole2.normalize();    
283 +
284 +    RealType cosAngle1 = dot(r12, r13);
285 +    RealType cosAngle2 = dot(r13, dipole2);
286 +
287 +    RealType halfBin = (nAngleBins_ - 1) * 0.5;
288 +    int angleBin1 = int(halfBin * (cosAngle1 + 1.0));
289 +    int angleBin2 = int(halfBin * (cosAngle2 + 1.0));
290 +
291 +    ++histogram_[angleBin1][angleBin2];    
292 +    ++npairs_;
293 +
294 +  }
295 +
296    void GofAngle2::writeRdf() {
297      std::ofstream rdfStream(outputFilename_.c_str());
298      if (rdfStream.is_open()) {
299        rdfStream << "#radial distribution function\n";
300        rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
301 <      rdfStream << "selection2: (" << selectionScript2_ << ")\n";
302 <      rdfStream << "#nAngleBins =" << nAngleBins_ << "deltaCosAngle = " << deltaCosAngle_ << "\n";
301 >      rdfStream << "selection2: (" << selectionScript2_ << ")";
302 >      if (doSele3_) {
303 >        rdfStream << "\tselection3: (" << selectionScript3_ << ")\n";
304 >      } else {
305 >        rdfStream << "\n";
306 >      }
307 >      rdfStream << "#nAngleBins =" << nAngleBins_ << "deltaCosAngle = "
308 >                << deltaCosAngle_ << "\n";
309        for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
310          // RealType cosAngle1 = -1.0 + (i + 0.5)*deltaCosAngle_;
311          
# Line 129 | Line 313 | namespace OpenMD {
313            // RealType cosAngle2 = -1.0 + (j + 0.5)*deltaCosAngle_;
314            rdfStream <<avgGofr_[i][j]/nProcessed_ << "\t";
315          }
132 –
316          rdfStream << "\n";
317        }
318          
319      } else {
320  
321 <      sprintf(painCave.errMsg, "GofAngle2: unable to open %s\n", outputFilename_.c_str());
321 >      sprintf(painCave.errMsg, "GofAngle2: unable to open %s\n",
322 >              outputFilename_.c_str());
323        painCave.isFatal = 1;
324        simError();  
325      }

Diff Legend

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