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

Comparing trunk/src/applications/staticProps/TetrahedralityParamZ.cpp (file contents):
Revision 1843 by gezelter, Tue Jan 29 20:58:08 2013 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   * [6]  Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
# Line 54 | Line 54 | namespace OpenMD {
54   namespace OpenMD {
55    TetrahedralityParamZ::TetrahedralityParamZ(SimInfo* info,  
56                                               const std::string& filename,
57 <                                             const std::string& sele,
57 >                                             const std::string& sele1,
58 >                                             const std::string& sele2,
59                                               double rCut, int nzbins)
60 <    : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info),
61 <      seleMan_(info), nZBins_(nzbins) {
60 >    : StaticAnalyser(info, filename),
61 >      selectionScript1_(sele1), selectionScript2_(sele2),
62 >      seleMan1_(info), seleMan2_(info), evaluator1_(info), evaluator2_(info),
63 >      nZBins_(nzbins) {
64      
65 <    evaluator_.loadScriptString(sele);
66 <    if (!evaluator_.isDynamic()) {
67 <      seleMan_.setSelectionSet(evaluator_.evaluate());
65 >    evaluator1_.loadScriptString(sele1);
66 >    if (!evaluator1_.isDynamic()) {
67 >      seleMan1_.setSelectionSet(evaluator1_.evaluate());
68      }
69 +    evaluator2_.loadScriptString(sele2);
70 +    if (!evaluator2_.isDynamic()) {
71 +      seleMan2_.setSelectionSet(evaluator2_.evaluate());
72 +    }
73      
74      // Set up cutoff radius:    
75      rCut_ = rCut;
# Line 91 | Line 98 | namespace OpenMD {
98      RigidBody* rb;
99      int myIndex;
100      SimInfo::MoleculeIterator mi;
94 –    Molecule::IntegrableObjectIterator ioi;
101      Molecule::RigidBodyIterator rbIter;
102      Vector3d vec;
103      Vector3d ri, rj, rk, rik, rkj;
# Line 99 | Line 105 | namespace OpenMD {
105      RealType cospsi;
106      RealType Qk;
107      std::vector<std::pair<RealType,StuntDouble*> > myNeighbors;
108 <    int isd;
108 >    int isd1;
109 >    int isd2;
110  
111      DumpReader reader(info_, dumpFilename_);    
112      int nFrames = reader.getNFrames();
# Line 113 | Line 120 | namespace OpenMD {
120        
121        RealType halfBoxZ_ = hmat(2,2) / 2.0;      
122  
123 <      if (evaluator_.isDynamic()) {
124 <        seleMan_.setSelectionSet(evaluator_.evaluate());
123 >      if (evaluator1_.isDynamic()) {
124 >        seleMan1_.setSelectionSet(evaluator1_.evaluate());
125        }
126        
127 +      if (evaluator2_.isDynamic()) {
128 +        seleMan2_.setSelectionSet(evaluator2_.evaluate());
129 +      }
130 +      
131        // update the positions of atoms which belong to the rigidbodies
132        for (mol = info_->beginMolecule(mi); mol != NULL;
133             mol = info_->nextMolecule(mi)) {
# Line 127 | Line 138 | namespace OpenMD {
138        }
139        
140        // outer loop is over the selected StuntDoubles:
141 <      for (sd = seleMan_.beginSelected(isd); sd != NULL;
142 <           sd = seleMan_.nextSelected(isd)) {
141 >      for (sd = seleMan1_.beginSelected(isd1); sd != NULL;
142 >           sd = seleMan1_.nextSelected(isd1)) {
143          
144          myIndex = sd->getGlobalIndex();
134 –
135 –        Qk = 1.0;        
136 –        myNeighbors.clear();
145          
146 <        // inner loop is over all StuntDoubles in the system:
147 <        
140 <        for (mol = info_->beginMolecule(mi); mol != NULL;
141 <             mol = info_->nextMolecule(mi)) {
146 >        Qk = 1.0;        
147 >        myNeighbors.clear();      
148  
149 <          for (sd2 = mol->beginIntegrableObject(ioi); sd2 != NULL;
150 <               sd2 = mol->nextIntegrableObject(ioi)) {
149 >        for (sd2 = seleMan2_.beginSelected(isd2); sd2 != NULL;
150 >             sd2 = seleMan2_.nextSelected(isd2)) {
151 >          
152 >          if (sd2->getGlobalIndex() != myIndex) {
153              
154 <            if (sd2->getGlobalIndex() != myIndex) {
155 <              
156 <              vec = sd->getPos() - sd2->getPos();      
157 <              
158 <              if (usePeriodicBoundaryConditions_)
159 <                currentSnapshot_->wrapVector(vec);
160 <              
161 <              r = vec.length();            
162 <
163 <              // Check to see if neighbor is in bond cutoff
164 <              
157 <              if (r < rCut_) {                
158 <                myNeighbors.push_back(std::make_pair(r,sd2));
159 <              }
154 >            vec = sd->getPos() - sd2->getPos();      
155 >            
156 >            if (usePeriodicBoundaryConditions_)
157 >              currentSnapshot_->wrapVector(vec);
158 >            
159 >            r = vec.length();            
160 >            
161 >            // Check to see if neighbor is in bond cutoff
162 >            
163 >            if (r < rCut_) {                
164 >              myNeighbors.push_back(std::make_pair(r,sd2));
165              }
166            }
167          }
# Line 170 | Line 175 | namespace OpenMD {
175          int nang = int (0.5 * (nbors * (nbors - 1)));
176          
177          rk = sd->getPos();
178 <
178 >        
179          for (int i = 0; i < nbors-1; i++) {      
180 <
180 >          
181            sdi = myNeighbors[i].second;
182            ri = sdi->getPos();
183            rik = rk - ri;
# Line 180 | Line 185 | namespace OpenMD {
185              currentSnapshot_->wrapVector(rik);
186            
187            rik.normalize();
188 <
188 >          
189            for (int j = i+1; j < nbors; j++) {      
190 <
190 >            
191              sdj = myNeighbors[j].second;
192              rj = sdj->getPos();
193              rkj = rk - rj;
# Line 191 | Line 196 | namespace OpenMD {
196              rkj.normalize();
197              
198              cospsi = dot(rik,rkj);          
199 <
199 >            
200              // Calculates scaled Qk for each molecule using calculated
201              // angles from 4 or fewer nearest neighbors.
202              Qk -=  (pow(cospsi + 1.0 / 3.0, 2) * 2.25 / nang);            
203            }
204          }
205 <
205 >        
206          if (nang > 0) {
207            if (usePeriodicBoundaryConditions_)
208              currentSnapshot_->wrapVector(rk);
# Line 226 | Line 231 | namespace OpenMD {
231      if (qZstream.is_open()) {
232        qZstream << "#Tetrahedrality Parameters (z)\n";
233        qZstream << "#nFrames:\t" << zBox_.size() << "\n";
234 <      qZstream << "#selection: (" << selectionScript_ << ")\n";
234 >      qZstream << "#selection 1: (" << selectionScript1_ << ")\n";
235 >      qZstream << "#selection 2: (" << selectionScript2_ << ")\n";
236        qZstream << "#z\tQk\n";
237        for (unsigned int i = 0; i < sliceQ_.size(); ++i) {
238          RealType z = zAve * (i+0.5) / sliceQ_.size();
239 <        qZstream << z << "\t" << sliceQ_[i] / sliceCount_[i] << "\n";
239 >        if (sliceCount_[i] != 0) {
240 >          qZstream << z << "\t" << sliceQ_[i] / sliceCount_[i] << "\n";
241 >        }
242        }
243        
244      } else {      

Diff Legend

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