ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/staticProps/TetrahedralityParam.cpp
Revision: 1785
Committed: Wed Aug 22 18:43:27 2012 UTC (12 years, 8 months ago) by jmichalk
File size: 10539 byte(s)
Log Message:
Trunk: The changes in this commit are confined to applications/staticProps and for the most part deal with a misspelling of initialize.

The one other change took place in StaticProps.cpp and deals with the default treatment of sele2. It had previously been set to 'select all' which seems to go against what would be desired by not specifying it with regard to proper operations of many of the analysis programs ( g of r's especially)

File Contents

# User Rev Content
1 kstocke1 1524 /*
2     * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3     *
4     * The University of Notre Dame grants you ("Licensee") a
5     * non-exclusive, royalty free, license to use, modify and
6     * redistribute this software in source and binary code form, provided
7     * that the following conditions are met:
8     *
9     * 1. Redistributions of source code must retain the above copyright
10     * notice, this list of conditions and the following disclaimer.
11     *
12     * 2. Redistributions in binary form must reproduce the above copyright
13     * notice, this list of conditions and the following disclaimer in the
14     * documentation and/or other materials provided with the
15     * distribution.
16     *
17     * This software is provided "AS IS," without a warranty of any
18     * kind. All express or implied conditions, representations and
19     * warranties, including any implied warranty of merchantability,
20     * fitness for a particular purpose or non-infringement, are hereby
21     * excluded. The University of Notre Dame and its licensors shall not
22     * be liable for any damages suffered by licensee as a result of
23     * using, modifying or distributing the software or its
24     * derivatives. In no event will the University of Notre Dame or its
25     * licensors be liable for any lost revenue, profit or data, or for
26     * direct, indirect, special, consequential, incidental or punitive
27     * damages, however caused and regardless of the theory of liability,
28     * arising out of the use of or inability to use software, even if the
29     * University of Notre Dame has been advised of the possibility of
30     * such damages.
31     *
32     * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33     * research, please cite the appropriate papers when you publish your
34     * work. Good starting points are:
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).
39 gezelter 1782 * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40     * [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
41 kstocke1 1524 * Created by J. Daniel Gezelter on 09/26/06.
42     * @author J. Daniel Gezelter
43     * @version $Id: BondOrderParameter.cpp 1442 2010-05-10 17:28:26Z gezelter $
44     *
45     */
46    
47     #include "applications/staticProps/TetrahedralityParam.hpp"
48     #include "utils/simError.h"
49     #include "io/DumpReader.hpp"
50     #include "primitives/Molecule.hpp"
51     #include "utils/NumericConstant.hpp"
52     #include <vector>
53    
54     namespace OpenMD {
55    
56     TetrahedralityParam::TetrahedralityParam(SimInfo* info,
57     const std::string& filename,
58     const std::string& sele,
59     double rCut, int nbins) : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info){
60    
61     setOutputName(getPrefix(filename) + ".q");
62    
63     evaluator_.loadScriptString(sele);
64     if (!evaluator_.isDynamic()) {
65     seleMan_.setSelectionSet(evaluator_.evaluate());
66     }
67    
68     // Set up cutoff radius:
69    
70     rCut_ = rCut;
71     nBins_ = nbins;
72    
73     Q_histogram_.resize(nBins_);
74    
75     // Q can take values from 0 to 1
76    
77     MinQ_ = 0.0;
78     MaxQ_ = 1.1;
79     deltaQ_ = (MaxQ_ - MinQ_) / nbins;
80    
81     }
82    
83     TetrahedralityParam::~TetrahedralityParam() {
84     Q_histogram_.clear();
85     }
86    
87 jmichalk 1785 void TetrahedralityParam::initializeHistogram() {
88 kstocke1 1524 std::fill(Q_histogram_.begin(), Q_histogram_.end(), 0);
89     }
90    
91     void TetrahedralityParam::process() {
92     Molecule* mol;
93     StuntDouble* sd;
94     StuntDouble* sd2;
95     StuntDouble* sdi;
96     StuntDouble* sdj;
97     RigidBody* rb;
98     int myIndex;
99     SimInfo::MoleculeIterator mi;
100     Molecule::RigidBodyIterator rbIter;
101     Molecule::IntegrableObjectIterator ioi;
102     Vector3d vec;
103     Vector3d ri, rj, rk, rik, rkj, dposition, tposition;
104     RealType r;
105     RealType cospsi;
106     RealType Qk;
107     std::vector<std::pair<RealType,StuntDouble*> > myNeighbors;
108     int isd;
109    
110     DumpReader reader(info_, dumpFilename_);
111     int nFrames = reader.getNFrames();
112     frameCounter_ = 0;
113    
114     Distorted_.clear();
115     Tetrahedral_.clear();
116    
117     for (int istep = 0; istep < nFrames; istep += step_) {
118     reader.readFrame(istep);
119     frameCounter_++;
120     currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
121    
122     if (evaluator_.isDynamic()) {
123     seleMan_.setSelectionSet(evaluator_.evaluate());
124     }
125    
126     // update the positions of atoms which belong to the rigidbodies
127    
128     for (mol = info_->beginMolecule(mi); mol != NULL;
129     mol = info_->nextMolecule(mi)) {
130     for (rb = mol->beginRigidBody(rbIter); rb != NULL;
131     rb = mol->nextRigidBody(rbIter)) {
132     rb->updateAtoms();
133     }
134     }
135    
136    
137     // outer loop is over the selected StuntDoubles:
138    
139     for (sd = seleMan_.beginSelected(isd); sd != NULL;
140     sd = seleMan_.nextSelected(isd)) {
141    
142     myIndex = sd->getGlobalIndex();
143     Qk = 1.0;
144    
145     myNeighbors.clear();
146    
147     // inner loop is over all StuntDoubles in the system:
148    
149     for (mol = info_->beginMolecule(mi); mol != NULL;
150     mol = info_->nextMolecule(mi)) {
151    
152     for (sd2 = mol->beginIntegrableObject(ioi); sd2 != NULL;
153     sd2 = mol->nextIntegrableObject(ioi)) {
154    
155     if (sd2->getGlobalIndex() != myIndex) {
156    
157     vec = sd->getPos() - sd2->getPos();
158    
159     if (usePeriodicBoundaryConditions_)
160     currentSnapshot_->wrapVector(vec);
161    
162     r = vec.length();
163    
164     // Check to see if neighbor is in bond cutoff
165    
166     if (r < rCut_) {
167    
168     myNeighbors.push_back(std::make_pair(r,sd2));
169     }
170     }
171     }
172     }
173    
174     // Sort the vector using predicate and std::sort
175     std::sort(myNeighbors.begin(), myNeighbors.end());
176    
177     std::cerr << myNeighbors.size() << " neighbors within " << rCut_ << " A" << " \n";
178    
179     // Use only the 4 closest neighbors to do the rest of the work:
180    
181 gezelter 1782 int nbors = myNeighbors.size()> 4 ? 4 : myNeighbors.size();
182 kstocke1 1524 int nang = int (0.5 * (nbors * (nbors - 1)));
183    
184     rk = sd->getPos();
185 gezelter 1782 std::cerr<<nbors<<endl;
186 kstocke1 1524 for (int i = 0; i < nbors-1; i++) {
187    
188     sdi = myNeighbors[i].second;
189     ri = sdi->getPos();
190     rik = rk - ri;
191     if (usePeriodicBoundaryConditions_)
192     currentSnapshot_->wrapVector(rik);
193    
194     rik.normalize();
195    
196     for (int j = i+1; j < nbors; j++) {
197    
198     sdj = myNeighbors[j].second;
199     rj = sdj->getPos();
200     rkj = rk - rj;
201     if (usePeriodicBoundaryConditions_)
202     currentSnapshot_->wrapVector(rkj);
203     rkj.normalize();
204    
205     cospsi = dot(rik,rkj);
206    
207     //std::cerr << "cos(psi) = " << cospsi << " \n";
208    
209     // Calculates scaled Qk for each molecule using calculated angles from 4 or fewer nearest neighbors.
210     Qk = Qk - (pow(cospsi + 1.0 / 3.0, 2) * 2.25 / nang);
211 gezelter 1782 //std::cerr<<Qk<<"\t"<<nang<<endl;
212 kstocke1 1524 }
213     }
214 gezelter 1782 //std::cerr<<nang<<endl;
215 kstocke1 1524 if (nang > 0) {
216     collectHistogram(Qk);
217    
218     // Saves positions of StuntDoubles & neighbors with distorted coordination (low Qk value)
219     if ((Qk < 0.55) && (Qk > 0.45)) {
220 gezelter 1782 //std::cerr<<Distorted_.size()<<endl;
221 kstocke1 1524 Distorted_.push_back(sd);
222 gezelter 1782 //std::cerr<<Distorted_.size()<<endl;
223 kstocke1 1524 dposition = sd->getPos();
224     //std::cerr << "distorted position \t" << dposition << "\n";
225     }
226    
227     // Saves positions of StuntDoubles & neighbors with tetrahedral coordination (high Qk value)
228 gezelter 1782 if (Qk > 0.05) {
229 kstocke1 1524
230     Tetrahedral_.push_back(sd);
231    
232     tposition = sd->getPos();
233     //std::cerr << "tetrahedral position \t" << tposition << "\n";
234     }
235    
236 gezelter 1782 //std::cerr<<Tetrahedral_.size()<<endl;
237    
238    
239 kstocke1 1524 }
240    
241     }
242     }
243    
244     writeOrderParameter();
245     std::cerr << "number of distorted StuntDoubles = " << Distorted_.size() << "\n";
246     std::cerr << "number of tetrahedral StuntDoubles = " << Tetrahedral_.size() << "\n";
247     }
248    
249     void TetrahedralityParam::collectHistogram(RealType Qk) {
250    
251     if (Qk > MinQ_ && Qk < MaxQ_) {
252    
253     int whichBin = int((Qk - MinQ_) / deltaQ_);
254     Q_histogram_[whichBin] += 1;
255     }
256     }
257    
258    
259     void TetrahedralityParam::writeOrderParameter() {
260    
261     int nSelected = 0;
262    
263     for (int i = 0; i < nBins_; ++i) {
264     nSelected = nSelected + Q_histogram_[i]*deltaQ_;
265     }
266 gezelter 1782
267 kstocke1 1524 std::ofstream osq((getOutputFileName() + "Q").c_str());
268    
269     if (osq.is_open()) {
270    
271     osq << "# Tetrahedrality Parameters\n";
272     osq << "# selection: (" << selectionScript_ << ")\n";
273     osq << "# \n";
274     // Normalize by number of frames and write it out:
275     for (int i = 0; i < nBins_; ++i) {
276     RealType Qval = MinQ_ + (i + 0.5) * deltaQ_;
277     osq << Qval;
278     osq << "\t" << (RealType) (Q_histogram_[i]/deltaQ_)/nSelected;
279     osq << "\n";
280     }
281    
282     osq.close();
283    
284     }else {
285     sprintf(painCave.errMsg, "TetrahedralityParam: unable to open %s\n",
286     (getOutputFileName() + "q").c_str());
287     painCave.isFatal = 1;
288     simError();
289     }
290    
291     DumpReader reader(info_, dumpFilename_);
292     int nFrames = reader.getNFrames();
293    
294     if (nFrames == 1) {
295    
296     std::vector<StuntDouble*>::iterator iter;
297     std::ofstream osd((getOutputFileName() + "dxyz").c_str());
298    
299     if (osd.is_open()) {
300    
301     osd << Distorted_.size() << "\n";
302    
303     osd << "1000000.00000000; 34.52893134 0.00000000 0.00000000; 0.00000000 34.52893134 0.00000000; 0.00000000 0.00000000 34.52893134" << "\n";
304    
305     for (iter = Distorted_.begin(); iter != Distorted_.end(); ++iter) {
306     Vector3d position;
307     position = (*iter)->getPos();
308     osd << "O " << "\t";
309 gezelter 1782 for (unsigned int z = 0; z < position.size(); z++) {
310 kstocke1 1524 osd << position[z] << " " << "\t";
311     }
312     osd << "\n";
313     }
314     osd.close();
315     }
316    
317    
318     std::ofstream ost((getOutputFileName() + "txyz").c_str());
319    
320     if (ost.is_open()) {
321    
322     ost << Tetrahedral_.size() << "\n";
323    
324     ost << "1000000.00000000; 34.52893134 0.00000000 0.00000000; 0.00000000 34.52893134 0.00000000; 0.00000000 0.00000000 34.52893134" << "\n";
325    
326     for (iter = Tetrahedral_.begin(); iter != Tetrahedral_.end(); ++iter) {
327    
328     Vector3d position;
329    
330     position = (*iter)->getPos();
331    
332     ost << "O " << "\t";
333    
334 gezelter 1782 for (unsigned int z = 0; z < position.size(); z++) {
335 kstocke1 1524
336     ost << position[z] << " " << "\t";
337     }
338    
339     ost << "\n";
340    
341     }
342    
343     ost.close();
344     }
345    
346     }
347     }
348     }
349    
350    
351    

Properties

Name Value
svn:executable *