ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/staticProps/TetrahedralityParam.cpp
Revision: 2071
Committed: Sat Mar 7 21:41:51 2015 UTC (10 years, 2 months ago) by gezelter
File size: 10290 byte(s)
Log Message:
Reducing the number of warnings when using g++ to compile.

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 gezelter 1879 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (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 gezelter 2071 double rCut, int nbins) :
60     StaticAnalyser(info, filename), selectionScript_(sele),
61     seleMan_(info), evaluator_(info) {
62 kstocke1 1524
63     setOutputName(getPrefix(filename) + ".q");
64 gezelter 2071
65 kstocke1 1524 evaluator_.loadScriptString(sele);
66     if (!evaluator_.isDynamic()) {
67     seleMan_.setSelectionSet(evaluator_.evaluate());
68     }
69    
70     // Set up cutoff radius:
71    
72     rCut_ = rCut;
73     nBins_ = nbins;
74    
75     Q_histogram_.resize(nBins_);
76    
77     // Q can take values from 0 to 1
78    
79     MinQ_ = 0.0;
80     MaxQ_ = 1.1;
81     deltaQ_ = (MaxQ_ - MinQ_) / nbins;
82    
83     }
84    
85     TetrahedralityParam::~TetrahedralityParam() {
86     Q_histogram_.clear();
87     }
88    
89 jmichalk 1785 void TetrahedralityParam::initializeHistogram() {
90 kstocke1 1524 std::fill(Q_histogram_.begin(), Q_histogram_.end(), 0);
91     }
92    
93     void TetrahedralityParam::process() {
94     Molecule* mol;
95     StuntDouble* sd;
96     StuntDouble* sd2;
97     StuntDouble* sdi;
98     StuntDouble* sdj;
99     RigidBody* rb;
100     int myIndex;
101     SimInfo::MoleculeIterator mi;
102     Molecule::RigidBodyIterator rbIter;
103     Molecule::IntegrableObjectIterator ioi;
104     Vector3d vec;
105     Vector3d ri, rj, rk, rik, rkj, dposition, tposition;
106     RealType r;
107     RealType cospsi;
108     RealType Qk;
109     std::vector<std::pair<RealType,StuntDouble*> > myNeighbors;
110     int isd;
111    
112     DumpReader reader(info_, dumpFilename_);
113     int nFrames = reader.getNFrames();
114     frameCounter_ = 0;
115    
116     Distorted_.clear();
117     Tetrahedral_.clear();
118    
119     for (int istep = 0; istep < nFrames; istep += step_) {
120     reader.readFrame(istep);
121     frameCounter_++;
122     currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
123    
124     if (evaluator_.isDynamic()) {
125     seleMan_.setSelectionSet(evaluator_.evaluate());
126     }
127    
128     // update the positions of atoms which belong to the rigidbodies
129    
130     for (mol = info_->beginMolecule(mi); mol != NULL;
131     mol = info_->nextMolecule(mi)) {
132     for (rb = mol->beginRigidBody(rbIter); rb != NULL;
133     rb = mol->nextRigidBody(rbIter)) {
134     rb->updateAtoms();
135     }
136     }
137    
138    
139 gezelter 2071 // outer loop is over the selected StuntDoubles:
140 kstocke1 1524
141     for (sd = seleMan_.beginSelected(isd); sd != NULL;
142     sd = seleMan_.nextSelected(isd)) {
143    
144     myIndex = sd->getGlobalIndex();
145     Qk = 1.0;
146    
147     myNeighbors.clear();
148    
149     // inner loop is over all StuntDoubles in the system:
150    
151     for (mol = info_->beginMolecule(mi); mol != NULL;
152     mol = info_->nextMolecule(mi)) {
153    
154     for (sd2 = mol->beginIntegrableObject(ioi); sd2 != NULL;
155     sd2 = mol->nextIntegrableObject(ioi)) {
156    
157     if (sd2->getGlobalIndex() != myIndex) {
158    
159     vec = sd->getPos() - sd2->getPos();
160    
161     if (usePeriodicBoundaryConditions_)
162     currentSnapshot_->wrapVector(vec);
163    
164     r = vec.length();
165    
166     // Check to see if neighbor is in bond cutoff
167    
168     if (r < rCut_) {
169    
170     myNeighbors.push_back(std::make_pair(r,sd2));
171     }
172     }
173     }
174     }
175    
176     // Sort the vector using predicate and std::sort
177     std::sort(myNeighbors.begin(), myNeighbors.end());
178    
179 gezelter 2071 //std::cerr << myNeighbors.size() << " neighbors within "
180     // << rCut_ << " A" << " \n";
181 kstocke1 1524
182     // Use only the 4 closest neighbors to do the rest of the work:
183    
184 gezelter 1782 int nbors = myNeighbors.size()> 4 ? 4 : myNeighbors.size();
185 kstocke1 1524 int nang = int (0.5 * (nbors * (nbors - 1)));
186    
187     rk = sd->getPos();
188 plouden 1786 //std::cerr<<nbors<<endl;
189 kstocke1 1524 for (int i = 0; i < nbors-1; i++) {
190    
191     sdi = myNeighbors[i].second;
192     ri = sdi->getPos();
193     rik = rk - ri;
194     if (usePeriodicBoundaryConditions_)
195     currentSnapshot_->wrapVector(rik);
196    
197     rik.normalize();
198    
199     for (int j = i+1; j < nbors; j++) {
200    
201     sdj = myNeighbors[j].second;
202     rj = sdj->getPos();
203     rkj = rk - rj;
204     if (usePeriodicBoundaryConditions_)
205     currentSnapshot_->wrapVector(rkj);
206     rkj.normalize();
207    
208     cospsi = dot(rik,rkj);
209    
210     //std::cerr << "cos(psi) = " << cospsi << " \n";
211    
212 gezelter 2071 // Calculates scaled Qk for each molecule using calculated
213     // angles from 4 or fewer nearest neighbors.
214 kstocke1 1524 Qk = Qk - (pow(cospsi + 1.0 / 3.0, 2) * 2.25 / nang);
215 gezelter 1782 //std::cerr<<Qk<<"\t"<<nang<<endl;
216 kstocke1 1524 }
217     }
218 gezelter 1782 //std::cerr<<nang<<endl;
219 kstocke1 1524 if (nang > 0) {
220     collectHistogram(Qk);
221    
222 gezelter 2071 // Saves positions of StuntDoubles & neighbors with distorted
223     // coordination (low Qk value)
224     if ((Qk < 0.55) && (Qk > 0.45)) {
225     //std::cerr<<Distorted_.size()<<endl;
226     Distorted_.push_back(sd);
227     //std::cerr<<Distorted_.size()<<endl;
228     dposition = sd->getPos();
229     //std::cerr << "distorted position \t" << dposition << "\n";
230     }
231 kstocke1 1524
232 gezelter 2071 // Saves positions of StuntDoubles & neighbors with
233     // tetrahedral coordination (high Qk value)
234     if (Qk > 0.05) {
235 kstocke1 1524
236 gezelter 2071 Tetrahedral_.push_back(sd);
237 kstocke1 1524
238 gezelter 2071 tposition = sd->getPos();
239     //std::cerr << "tetrahedral position \t" << tposition << "\n";
240     }
241 kstocke1 1524
242 gezelter 2071 //std::cerr<<Tetrahedral_.size()<<endl;
243    
244 kstocke1 1524 }
245    
246     }
247     }
248    
249     writeOrderParameter();
250 gezelter 2071 std::cerr << "number of distorted StuntDoubles = "
251     << Distorted_.size() << "\n";
252     std::cerr << "number of tetrahedral StuntDoubles = "
253     << Tetrahedral_.size() << "\n";
254 kstocke1 1524 }
255    
256     void TetrahedralityParam::collectHistogram(RealType Qk) {
257    
258     if (Qk > MinQ_ && Qk < MaxQ_) {
259    
260     int whichBin = int((Qk - MinQ_) / deltaQ_);
261     Q_histogram_[whichBin] += 1;
262     }
263     }
264    
265    
266     void TetrahedralityParam::writeOrderParameter() {
267    
268     int nSelected = 0;
269    
270     for (int i = 0; i < nBins_; ++i) {
271 gezelter 1790 nSelected = nSelected + int(Q_histogram_[i]*deltaQ_);
272 kstocke1 1524 }
273 gezelter 1782
274 kstocke1 1524 std::ofstream osq((getOutputFileName() + "Q").c_str());
275    
276     if (osq.is_open()) {
277    
278     osq << "# Tetrahedrality Parameters\n";
279     osq << "# selection: (" << selectionScript_ << ")\n";
280     osq << "# \n";
281     // Normalize by number of frames and write it out:
282     for (int i = 0; i < nBins_; ++i) {
283     RealType Qval = MinQ_ + (i + 0.5) * deltaQ_;
284     osq << Qval;
285     osq << "\t" << (RealType) (Q_histogram_[i]/deltaQ_)/nSelected;
286     osq << "\n";
287     }
288    
289     osq.close();
290    
291     }else {
292     sprintf(painCave.errMsg, "TetrahedralityParam: unable to open %s\n",
293     (getOutputFileName() + "q").c_str());
294     painCave.isFatal = 1;
295     simError();
296     }
297    
298     DumpReader reader(info_, dumpFilename_);
299     int nFrames = reader.getNFrames();
300    
301     if (nFrames == 1) {
302    
303 gezelter 2071 std::vector<StuntDouble*>::iterator iter;
304     std::ofstream osd((getOutputFileName() + "dxyz").c_str());
305 kstocke1 1524
306 gezelter 2071 if (osd.is_open()) {
307 kstocke1 1524
308 gezelter 2071 osd << Distorted_.size() << "\n";
309     osd << "\n";
310 kstocke1 1524
311 gezelter 2071 for (iter = Distorted_.begin(); iter != Distorted_.end(); ++iter) {
312     Vector3d position;
313     position = (*iter)->getPos();
314     osd << "O " << "\t";
315 gezelter 1782 for (unsigned int z = 0; z < position.size(); z++) {
316 kstocke1 1524 osd << position[z] << " " << "\t";
317     }
318     osd << "\n";
319 gezelter 2071 }
320     osd.close();
321 kstocke1 1524 }
322    
323    
324 gezelter 2071 std::ofstream ost((getOutputFileName() + "txyz").c_str());
325 kstocke1 1524
326 gezelter 2071 if (ost.is_open()) {
327 kstocke1 1524
328 gezelter 2071 ost << Tetrahedral_.size() << "\n";
329     ost << "\n";
330 kstocke1 1524
331 gezelter 2071 for (iter = Tetrahedral_.begin(); iter != Tetrahedral_.end(); ++iter) {
332     Vector3d position;
333     position = (*iter)->getPos();
334 kstocke1 1524
335 gezelter 2071 ost << "O " << "\t";
336 kstocke1 1524
337 gezelter 1782 for (unsigned int z = 0; z < position.size(); z++) {
338 kstocke1 1524 ost << position[z] << " " << "\t";
339     }
340     ost << "\n";
341 gezelter 2071 }
342     ost.close();
343 kstocke1 1524 }
344     }
345     }
346     }
347    
348    
349 gezelter 2071

Properties

Name Value
svn:executable *