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 1795 by gezelter, Wed Aug 22 02:28:28 2012 UTC vs.
Revision 1796 by gezelter, Mon Sep 10 18:38:44 2012 UTC

# Line 112 | Line 112 | namespace OpenMD
112      int myIndex;
113      SimInfo::MoleculeIterator mi;
114      Molecule::RigidBodyIterator rbIter;
115 –    Molecule::IntegrableObjectIterator ioi;
115      Vector3d vec;
116      Vector3d ri, rj, rk, rik, rkj, dposition, tposition;
117      RealType r;
# Line 138 | Line 137 | namespace OpenMD
137      Distorted_.clear();
138      Tetrahedral_.clear();
139      int i;
140 <    for(i=0;i<nZBins_;i++)
141 <      {
142 <        sliceSDLists_[i].clear();
144 <      }
140 >    for(i=0;i<nZBins_;i++) {
141 >      sliceSDLists_[i].clear();
142 >    }
143  
144      //LOOP OVER ALL FRAMES
145 <    for (int istep = 0; istep < nFrames; istep += step_)
146 <      {
147 <        int i;
148 <        for(i=0;i<nZBins_;i++)
149 <          {
152 <            count_[i]=0;
153 <          }    
145 >    for (int istep = 0; istep < nFrames; istep += step_) {
146 >      int i;
147 >      for(i=0;i<nZBins_;i++) {
148 >        count_[i]=0;
149 >      }
150          
151 <        reader.readFrame(istep);
152 <        frameCounter_++;
153 <        currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
151 >      reader.readFrame(istep);
152 >      frameCounter_++;
153 >      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
154 >      
155 >      if (evaluator_.isDynamic()) {
156 >        seleMan1_.setSelectionSet(evaluator_.evaluate());
157 >        seleMan2_.setSelectionSet(evaluator_.evaluate());
158 >      }
159          
159 –        if (evaluator_.isDynamic())
160 –          {
161 –            seleMan1_.setSelectionSet(evaluator_.evaluate());
162 –            seleMan2_.setSelectionSet(evaluator_.evaluate());
163 –          }
164 –        
160        // update the positions of atoms which belong to the rigidbodies
161 <      for (mol = info_->beginMolecule(mi); mol != NULL;mol = info_->nextMolecule(mi))
162 <        {
163 <          for (rb = mol->beginRigidBody(rbIter); rb != NULL;rb = mol->nextRigidBody(rbIter))
164 <            {
165 <              rb->updateAtoms();
166 <            }        
167 <        }          
161 >      for (mol = info_->beginMolecule(mi); mol != NULL;
162 >           mol = info_->nextMolecule(mi)) {
163 >        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
164 >             rb = mol->nextRigidBody(rbIter)) {
165 >          rb->updateAtoms();
166 >        }        
167 >      }          
168  
169 <       // outer loop is over the selected StuntDoubles:
169 >      // outer loop is over the selected StuntDoubles:
170        int idk=0;
171 <      for (sd = seleMan1_.beginSelected(isd1); sd != NULL;sd = seleMan1_.nextSelected(isd1))
172 <        {
173 <          myIndex = sd->getGlobalIndex();
174 <          Qk = 1.0;      
175 <          myNeighbors.clear();
176 <          // inner loop is over all StuntDoubles in the system:
177 <          //for (mol = info_->beginMolecule(mi); mol != NULL;mol = info_->nextMolecule(mi))
178 <            //{
179 <              //for (sd2 = mol->beginIntegrableObject(ioi); sd2 != NULL; sd2 = mol->nextIntegrableObject(ioi))
180 <                //{
181 <                for(sd2 = seleMan2_.beginSelected(isd2); sd2 != NULL; sd2 = seleMan2_.nextSelected(isd2)){
182 <                        if(sd2->getGlobalIndex() != myIndex){
183 <                        vec = sd->getPos() - sd2->getPos();      
184 <                        if (usePeriodicBoundaryConditions_)
185 <                                        currentSnapshot_->wrapVector(vec);
186 <                        r = vec.length();            
187 <
188 <                        // Check to see if neighbor is in bond cutoff
189 <                        if (r < rCut_){
190 <                                        myNeighbors.push_back(std::make_pair(r,sd2));
191 <                                }
192 <                        }
193 <                }
194 <           // }
195 <          // Sort the vector using predicate and std::sort
196 <          std::sort(myNeighbors.begin(), myNeighbors.end());      
197 <          //std::cerr << myNeighbors.size() <<  " neighbors within " << rCut_  << " A" << " \n";
198 <          // Use only the 4 closest neighbors to do the rest of the work:        
199 <          int nbors =  myNeighbors.size()> 4 ? 4 : myNeighbors.size();
200 <          int nang = int (0.5 * (nbors * (nbors - 1)));
201 <          
202 <          rk = sd->getPos();
203 <          for (int i = 0; i < nbors-1; i++)
204 <            {                
205 <              sdi = myNeighbors[i].second;
206 <              ri = sdi->getPos();
207 <              rik = rk - ri;
208 <              if (usePeriodicBoundaryConditions_)
209 <                currentSnapshot_->wrapVector(rik);            
210 <              rik.normalize();
211 <              
212 <              for (int j = i+1; j < nbors; j++)
213 <                {                    
214 <                  sdj = myNeighbors[j].second;
215 <                  rj = sdj->getPos();
216 <                  rkj = rk - rj;
217 <                  if (usePeriodicBoundaryConditions_)
218 <                    currentSnapshot_->wrapVector(rkj);
219 <                  rkj.normalize();
220 <                  
221 <                  cospsi = dot(rik,rkj);
222 <                  
223 <                  // Calculates scaled Qk for each molecule using calculated angles from 4 or fewer nearest neighbors.
224 <                  Qk = Qk - (pow(cospsi + 1.0 / 3.0, 2) * 2.25 / nang);
225 <                }
226 <            }
227 <          
228 <          //std::cerr<<nbors<<endl;
229 <           if (nang > 0)
230 <           {
231 <             //collectHistogram(Qk);
232 <              
233 <              // Saves positions of StuntDoubles & neighbors with distorted coordination (low Qk value)
234 <            if ((Qk < 0.55) && (Qk > 0.45))
235 <                {
236 <                  Distorted_.push_back(sd);
237 <                  dposition = sd->getPos();
238 <                }
239 <            
240 <            // Saves positions of StuntDoubles & neighbors with tetrahedral coordination (high Qk value)
241 <            if (Qk > 0)
242 <              {
243 <                Tetrahedral_.push_back(sd);
244 <                tposition = sd->getPos();
245 <              }
246 <
247 <           }
248 <          
249 <           //wrap the stuntdoubles into a cell      
250 <           Vector3d pos = sd->getPos();
251 <           if (usePeriodicBoundaryConditions_)
252 <             currentSnapshot_->wrapVector(pos);
253 <           sd->setPos(pos);
254 <           // shift molecules by half a box to have bins start at 0
255 <           int binNo = int(nZBins_ * (halfBoxZ_ + pos.z()) / hmat(2,2));
261 <           //Patrick took out the "halfBoxZ_" part in the line above to below
262 <           //int binNo = int(nZBins_ * (pos.z()) / hmat(2,2));
263 <           sliceSDLists_[binNo].push_back(Qk);
264 <           idk++;
265 <        }//outer sd loop
266 <      }//istep loop
171 >      for (sd = seleMan1_.beginSelected(isd1); sd != NULL;
172 >           sd = seleMan1_.nextSelected(isd1)) {
173 >        myIndex = sd->getGlobalIndex();
174 >        Qk = 1.0;        
175 >        myNeighbors.clear();
176 >        for(sd2 = seleMan2_.beginSelected(isd2); sd2 != NULL;
177 >            sd2 = seleMan2_.nextSelected(isd2)){
178 >          if(sd2->getGlobalIndex() != myIndex){
179 >            vec = sd->getPos() - sd2->getPos();      
180 >            if (usePeriodicBoundaryConditions_)
181 >              currentSnapshot_->wrapVector(vec);
182 >            r = vec.length();            
183 >            
184 >            // Check to see if neighbor is in bond cutoff
185 >            if (r < rCut_) {
186 >              myNeighbors.push_back(std::make_pair(r,sd2));
187 >            }
188 >          }
189 >        }
190 >        
191 >        // Sort the vector using predicate and std::sort
192 >        std::sort(myNeighbors.begin(), myNeighbors.end());        
193 >        //std::cerr << myNeighbors.size() <<  " neighbors within " << rCut_  << " A" << " \n";
194 >        // Use only the 4 closest neighbors to do the rest of the work:  
195 >        int nbors =  myNeighbors.size() > 4 ? 4 : myNeighbors.size();
196 >        int nang = int (0.5 * (nbors * (nbors - 1)));
197 >        
198 >        rk = sd->getPos();
199 >        for (int i = 0; i < nbors-1; i++) {
200 >          sdi = myNeighbors[i].second;
201 >          ri = sdi->getPos();
202 >          rik = rk - ri;
203 >          if (usePeriodicBoundaryConditions_)
204 >            currentSnapshot_->wrapVector(rik);        
205 >          rik.normalize();
206 >          
207 >          for (int j = i+1; j < nbors; j++) {
208 >            sdj = myNeighbors[j].second;
209 >            rj = sdj->getPos();
210 >            rkj = rk - rj;
211 >            if (usePeriodicBoundaryConditions_)
212 >              currentSnapshot_->wrapVector(rkj);
213 >            rkj.normalize();
214 >            
215 >            cospsi = dot(rik,rkj);
216 >            
217 >            // Calculates scaled Qk for each molecule using calculated
218 >            // angles from 4 or fewer nearest neighbors.
219 >            Qk = Qk - (pow(cospsi + 1.0 / 3.0, 2) * 2.25 / nang);
220 >          }
221 >        }
222 >        
223 >        //std::cerr<<nbors<<endl;
224 >        if (nang > 0) {
225 >          //collectHistogram(Qk);
226 >          
227 >          // Saves positions of StuntDoubles & neighbors with
228 >          // distorted coordination (low Qk value)
229 >          if ((Qk < 0.55) && (Qk > 0.45)) {
230 >            Distorted_.push_back(sd);
231 >            dposition = sd->getPos();
232 >          }
233 >          
234 >          // Saves positions of StuntDoubles & neighbors with
235 >          // tetrahedral coordination (high Qk value)
236 >          if (Qk > 0) {
237 >            Tetrahedral_.push_back(sd);
238 >            tposition = sd->getPos();
239 >          }
240 >          
241 >        }
242 >        
243 >        //wrap the stuntdoubles into a cell      
244 >        Vector3d pos = sd->getPos();
245 >        if (usePeriodicBoundaryConditions_)
246 >          currentSnapshot_->wrapVector(pos);
247 >        sd->setPos(pos);
248 >        // shift molecules by half a box to have bins start at 0
249 >        int binNo = int(nZBins_ * (halfBoxZ_ + pos.z()) / hmat(2,2));
250 >        // Patrick took out the "halfBoxZ_" part in the line above to below
251 >        // int binNo = int(nZBins_ * (pos.z()) / hmat(2,2));
252 >        sliceSDLists_[binNo].push_back(Qk);
253 >        idk++;
254 >      }//outer sd loop
255 >    }//istep loop
256  
257      //Averaging the value of Qk in each bin
258 <    for(int i=0;i< nZBins_; i++)
259 <      {
260 <        RealType Qsum=0;
261 <        for (unsigned int k = 0; k < sliceSDLists_[i].size(); ++k)
262 <          {
274 <            Qsum=Qsum+sliceSDLists_[i][k];
275 <            count_[i]++;
276 <          }
277 <        //std::cerr<<"past averagin Qk"<<endl;
278 <        //std::cerr<<Qsum<<endl;
279 <        if(count_[i]!=0)
280 <          {
281 <            Qave_.push_back(Qsum/count_[i]);
282 <          }
283 <        //std::cerr<<count[i]<<endl;
258 >    for(int i=0; i< nZBins_; i++) {
259 >      RealType Qsum=0;
260 >      for (unsigned int k = 0; k < sliceSDLists_[i].size(); ++k) {        
261 >        Qsum=Qsum+sliceSDLists_[i][k];
262 >        count_[i]++;
263        }
264 +      //std::cerr<<"past averagin Qk"<<endl;
265 +      //std::cerr<<Qsum<<endl;
266 +      if(count_[i]!=0) {
267 +        Qave_.push_back(Qsum/count_[i]);
268 +      }
269 +      //std::cerr<<count[i]<<endl;
270 +    }
271      //std::cerr<<"nZBins_ = "<< nZBins_<<endl;
272      //Writing bin#:<Qk> to a file
273      std::ofstream rdfStream(outputFilename_.c_str());
274 <    if (rdfStream.is_open())
275 <      {
276 <        //rdfStream << "#QkZ\n";
277 <        //rdfStream << "#nFrames:\t" << nProcessed_ << "\n";
278 <        //rdfStream << "#selection: (" << selectionScript_ << ")\n";
279 <        //rdfStream << "#z\tdensity\n";
280 <        for (int i = 0; i < nZBins_; ++i)
281 <          {
282 <            if(count_[i]!=0)
283 <              {
298 <                rdfStream << ((hmat(2,2)*i)/nZBins_)+(hmat(2,2)/(2*nZBins_)) << "\t" << Qave_[i] << "\n";
299 <              }
300 <          }
274 >    if (rdfStream.is_open()) {
275 >      //rdfStream << "#QkZ\n";
276 >      //rdfStream << "#nFrames:\t" << nProcessed_ << "\n";
277 >      //rdfStream << "#selection: (" << selectionScript_ << ")\n";
278 >      //rdfStream << "#z\tdensity\n";
279 >      for (int i = 0; i < nZBins_; ++i) {
280 >        if(count_[i]!=0) {
281 >          rdfStream << ((hmat(2,2)*i)/nZBins_)+(hmat(2,2)/(2*nZBins_))
282 >                    << "\t" << Qave_[i] << "\n";
283 >        }
284        }
285 +    }
286      
303 –
304 –
305 –
306 –
307 –  
287      writeOrderParameter();
288 <    std::cerr << "number of distorted StuntDoubles = " << Distorted_.size() << "\n";
289 <    std::cerr << "number of tetrahedral StuntDoubles = " << Tetrahedral_.size() << "\n";
288 >    std::cerr << "number of distorted StuntDoubles = "
289 >              << Distorted_.size() << "\n";
290 >    std::cerr << "number of tetrahedral StuntDoubles = "
291 >              << Tetrahedral_.size() << "\n";
292      collectHistogram(Qk);
293  
294    }//void TetrahedralityParam::process() loop
295    
296 <  void TetrahedralityParamZ::collectHistogram(RealType Qk)
297 <  {
298 <  //if (Qk > MinQ_ && Qk < MaxQ_)
299 <  //  {
300 <  //    int whichBin = int((Qk - MinQ_) / deltaQ_);
301 <  //    Q_histogram_[whichBin] += 1;
321 <  //  }
296 >  void TetrahedralityParamZ::collectHistogram(RealType Qk) {
297 >    //if (Qk > MinQ_ && Qk < MaxQ_)
298 >    //  {
299 >    //  int whichBin = int((Qk - MinQ_) / deltaQ_);
300 >    //  Q_histogram_[whichBin] += 1;
301 >    //  }
302    }    
303  
304 <  void TetrahedralityParamZ::writeOrderParameter()
305 <  {  
306 <   int nSelected = 0;
307 <  std::ofstream osq((getOutputFileName() + "Q").c_str());
308 <  if (osq.is_open())
309 <    {
310 <        osq << "# Tetrahedrality Parameters\n";
311 <        osq << "# selection: (" << selectionScript_ << ")\n";
312 <        osq << "# \n";
313 <        osq.close();
304 >  void TetrahedralityParamZ::writeOrderParameter() {
305 >    int nSelected = 0;
306 >    std::ofstream osq((getOutputFileName() + "Q").c_str());
307 >    if (osq.is_open()) {
308 >      osq << "# Tetrahedrality Parameters\n";
309 >      osq << "# selection: (" << selectionScript_ << ")\n";
310 >      osq << "# \n";
311 >      osq.close();
312 >    } else {
313 >      sprintf(painCave.errMsg, "TetrahedralityParamZ: unable to open %s\n",
314 >              (getOutputFileName() + "q").c_str());
315 >      painCave.isFatal = 1;
316 >      simError();  
317      }
318 <  else
319 <    {
320 <        sprintf(painCave.errMsg, "TetrahedralityParamZ: unable to open %s\n",
321 <                (getOutputFileName() + "q").c_str());
322 <        painCave.isFatal = 1;
323 <        simError();  
318 >    DumpReader reader(info_, dumpFilename_);    
319 >    int nFrames = reader.getNFrames();
320 >    if (nFrames == 1) {
321 >      std::vector<StuntDouble*>::iterator iter;
322 >      std::ofstream osd((getOutputFileName() + "dxyz").c_str());
323 >      if (osd.is_open()) {
324 >        osd << Distorted_.size() << "\n\n";
325 >        
326 >        for (iter = Distorted_.begin(); iter != Distorted_.end(); ++iter) {
327 >          Vector3d position;
328 >          position = (*iter)->getPos();
329 >          osd << "O  " << "\t";
330 >          for (unsigned int z=0; z<position.size(); z++) {
331 >            osd << position[z] << "  " << "\t";
332 >          }
333 >          osd << "\n";
334 >        }
335 >        osd.close();
336 >      }
337 >      std::ofstream ost((getOutputFileName() + "txyz").c_str());
338 >      if (ost.is_open()) {
339 >        ost << Tetrahedral_.size() << "\n\n";      
340 >        for (iter = Tetrahedral_.begin(); iter != Tetrahedral_.end(); ++iter) {
341 >          Vector3d position;            
342 >          position = (*iter)->getPos();
343 >          ost << "O  " << "\t";
344 >          for (unsigned int z=0; z<position.size(); z++) {
345 >            ost << position[z] << "  " << "\t";
346 >          }
347 >          ost << "\n";
348 >        }
349 >        ost.close();
350 >      }
351      }
342 –  DumpReader reader(info_, dumpFilename_);    
343 –  int nFrames = reader.getNFrames();
344 –  if (nFrames == 1)
345 –    {
346 –        std::vector<StuntDouble*>::iterator iter;
347 –        std::ofstream osd((getOutputFileName() + "dxyz").c_str());
348 –        if (osd.is_open())
349 –          {
350 –            osd << Distorted_.size() << "\n";
351 –            osd << "1000000.00000000;    34.52893134     0.00000000     0.00000000;     0.00000000    34.52893134     0.00000000;     0.00000000     0.00000000    34.52893134" << "\n";
352 –            
353 –            for (iter = Distorted_.begin(); iter != Distorted_.end(); ++iter)
354 –              {
355 –                Vector3d position;
356 –                position = (*iter)->getPos();
357 –                osd << "O  " << "\t";
358 –                for (unsigned int z=0; z<position.size(); z++)
359 –                  {
360 –                    osd << position[z] << "  " << "\t";
361 –                  }
362 –                osd << "\n";
363 –              }
364 –            osd.close();
365 –          }
366 –        std::ofstream ost((getOutputFileName() + "txyz").c_str());
367 –        if (ost.is_open())
368 –          {
369 –            ost << Tetrahedral_.size() << "\n";
370 –            ost << "1000000.00000000;    34.52893134     0.00000000     0.00000000;     0.00000000    34.52893134     0.00000000;     0.00000000     0.00000000    34.52893134" << "\n";
371 –            
372 –            for (iter = Tetrahedral_.begin(); iter != Tetrahedral_.end(); ++iter)
373 –              {
374 –                Vector3d position;              
375 –                position = (*iter)->getPos();
376 –                ost << "O  " << "\t";
377 –                for (unsigned int z=0; z<position.size(); z++)
378 –                  {
379 –                    ost << position[z] << "  " << "\t";
380 –                  }
381 –                ost << "\n";
382 –              }
383 –            ost.close();
384 –          }
385 –    }
352    }
353   }
354  

Diff Legend

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