OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
TetrahedralityParam.cpp
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the appropriate papers when you publish your
33 * work. Good starting points are:
34 *
35 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
38 * [4] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
39 * [5] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
40 * [6] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
41 * [7] Lamichhane, Newman & Gezelter, J. Chem. Phys. 141, 134110 (2014).
42 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
43 */
44
45#include "applications/staticProps/TetrahedralityParam.hpp"
46
47#include <vector>
48
49#include "io/DumpReader.hpp"
51#include "utils/simError.h"
52
53namespace OpenMD {
54
55 TetrahedralityParam::TetrahedralityParam(SimInfo* info,
56 const std::string& filename,
57 const std::string& sele, double rCut,
58 int nbins) :
59 StaticAnalyser(info, filename, nbins),
60 selectionScript_(sele), seleMan_(info), evaluator_(info) {
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
72 Q_histogram_.resize(nBins_);
73
74 // Q can take values from 0 to 1
75
76 MinQ_ = 0.0;
77 MaxQ_ = 1.1;
78 deltaQ_ = (MaxQ_ - MinQ_) / nBins_;
79 }
80
81 void TetrahedralityParam::initializeHistogram() {
82 std::fill(Q_histogram_.begin(), Q_histogram_.end(), 0);
83 }
84
85 void TetrahedralityParam::process() {
86 Molecule* mol;
87 StuntDouble* sd;
88 StuntDouble* sd2;
89 StuntDouble* sdi;
90 StuntDouble* sdj;
91 int myIndex;
92 SimInfo::MoleculeIterator mi;
93 Molecule::IntegrableObjectIterator ioi;
94 Vector3d vec;
95 Vector3d ri, rj, rk, rik, rkj, dposition, tposition;
96 RealType r;
97 RealType cospsi;
98 RealType Qk;
99 std::vector<std::pair<RealType, StuntDouble*>> myNeighbors;
100 int isd;
101 bool usePeriodicBoundaryConditions_ =
102 info_->getSimParams()->getUsePeriodicBoundaryConditions();
103
104 DumpReader reader(info_, dumpFilename_);
105 int nFrames = reader.getNFrames();
106 frameCounter_ = 0;
107
108 Distorted_.clear();
109 Tetrahedral_.clear();
110
111 for (int istep = 0; istep < nFrames; istep += step_) {
112 reader.readFrame(istep);
113 frameCounter_++;
114 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
115
116 if (evaluator_.isDynamic()) {
117 seleMan_.setSelectionSet(evaluator_.evaluate());
118 }
119
120 // outer loop is over the selected StuntDoubles:
121
122 for (sd = seleMan_.beginSelected(isd); sd != NULL;
123 sd = seleMan_.nextSelected(isd)) {
124 myIndex = sd->getGlobalIndex();
125 Qk = 1.0;
126
127 myNeighbors.clear();
128
129 // inner loop is over all StuntDoubles in the system:
130
131 for (mol = info_->beginMolecule(mi); mol != NULL;
132 mol = info_->nextMolecule(mi)) {
133 for (sd2 = mol->beginIntegrableObject(ioi); sd2 != NULL;
134 sd2 = mol->nextIntegrableObject(ioi)) {
135 if (sd2->getGlobalIndex() != myIndex) {
136 vec = sd->getPos() - sd2->getPos();
137
138 if (usePeriodicBoundaryConditions_)
139 currentSnapshot_->wrapVector(vec);
140
141 r = vec.length();
142
143 // Check to see if neighbor is in bond cutoff
144
145 if (r < rCut_) { myNeighbors.push_back(std::make_pair(r, sd2)); }
146 }
147 }
148 }
149
150 // Sort the vector using predicate and std::sort
151 std::sort(myNeighbors.begin(), myNeighbors.end());
152
153 // std::cerr << myNeighbors.size() << " neighbors within "
154 // << rCut_ << " A" << " \n";
155
156 // Use only the 4 closest neighbors to do the rest of the work:
157
158 int nbors = myNeighbors.size() > 4 ? 4 : myNeighbors.size();
159 int nang = int(0.5 * (nbors * (nbors - 1)));
160
161 rk = sd->getPos();
162 // std::cerr<<nbors<<endl;
163 for (int i = 0; i < nbors - 1; i++) {
164 sdi = myNeighbors[i].second;
165 ri = sdi->getPos();
166 rik = rk - ri;
167 if (usePeriodicBoundaryConditions_) currentSnapshot_->wrapVector(rik);
168
169 rik.normalize();
170
171 for (int j = i + 1; j < nbors; j++) {
172 sdj = myNeighbors[j].second;
173 rj = sdj->getPos();
174 rkj = rk - rj;
175 if (usePeriodicBoundaryConditions_)
176 currentSnapshot_->wrapVector(rkj);
177 rkj.normalize();
178
179 cospsi = dot(rik, rkj);
180
181 // std::cerr << "cos(psi) = " << cospsi << " \n";
182
183 // Calculates scaled Qk for each molecule using calculated
184 // angles from 4 or fewer nearest neighbors.
185 Qk = Qk - (pow(cospsi + 1.0 / 3.0, 2) * 2.25 / nang);
186 // std::cerr<<Qk<<"\t"<<nang<<endl;
187 }
188 }
189 // std::cerr<<nang<<endl;
190 if (nang > 0) {
191 collectHistogram(Qk);
192
193 // Saves positions of StuntDoubles & neighbors with distorted
194 // coordination (low Qk value)
195 if ((Qk < 0.55) && (Qk > 0.45)) {
196 // std::cerr<<Distorted_.size()<<endl;
197 Distorted_.push_back(sd);
198 // std::cerr<<Distorted_.size()<<endl;
199 dposition = sd->getPos();
200 // std::cerr << "distorted position \t" << dposition << "\n";
201 }
202
203 // Saves positions of StuntDoubles & neighbors with
204 // tetrahedral coordination (high Qk value)
205 if (Qk > 0.05) {
206 Tetrahedral_.push_back(sd);
207
208 tposition = sd->getPos();
209 // std::cerr << "tetrahedral position \t" << tposition << "\n";
210 }
211
212 // std::cerr<<Tetrahedral_.size()<<endl;
213 }
214 }
215 }
216
217 writeOrderParameter();
218 std::cerr << "number of distorted StuntDoubles = " << Distorted_.size()
219 << "\n";
220 std::cerr << "number of tetrahedral StuntDoubles = " << Tetrahedral_.size()
221 << "\n";
222 }
223
224 void TetrahedralityParam::collectHistogram(RealType Qk) {
225 if (Qk > MinQ_ && Qk < MaxQ_) {
226 int whichBin = int((Qk - MinQ_) / deltaQ_);
227 Q_histogram_[whichBin] += 1;
228 }
229 }
230
231 void TetrahedralityParam::writeOrderParameter() {
232 int nSelected = 0;
233
234 for (int i = 0; i < nBins_; ++i) {
235 nSelected = nSelected + int(Q_histogram_[i] * deltaQ_);
236 }
237
238 std::ofstream osq((getOutputFileName() + "Q").c_str());
239
240 if (osq.is_open()) {
241 osq << "# Tetrahedrality Parameters\n";
242 osq << "# selection: (" << selectionScript_ << ")\n";
243 osq << "# \n";
244 // Normalize by number of frames and write it out:
245 for (int i = 0; i < nBins_; ++i) {
246 RealType Qval = MinQ_ + (i + 0.5) * deltaQ_;
247 osq << Qval;
248 osq << "\t" << (RealType)(Q_histogram_[i] / deltaQ_) / nSelected;
249 osq << "\n";
250 }
251
252 osq.close();
253
254 } else {
255 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
256 "TetrahedralityParam: unable to open %s\n",
257 (getOutputFileName() + "q").c_str());
258 painCave.isFatal = 1;
259 simError();
260 }
261
262 DumpReader reader(info_, dumpFilename_);
263 int nFrames = reader.getNFrames();
264
265 if (nFrames == 1) {
266 std::vector<StuntDouble*>::iterator iter;
267 std::ofstream osd((getOutputFileName() + "dxyz").c_str());
268
269 if (osd.is_open()) {
270 osd << Distorted_.size() << "\n";
271 osd << "\n";
272
273 for (iter = Distorted_.begin(); iter != Distorted_.end(); ++iter) {
274 Vector3d position;
275 position = (*iter)->getPos();
276 osd << "O "
277 << "\t";
278 for (unsigned int z = 0; z < position.size(); z++) {
279 osd << position[z] << " "
280 << "\t";
281 }
282 osd << "\n";
283 }
284 osd.close();
285 }
286
287 std::ofstream ost((getOutputFileName() + "txyz").c_str());
288
289 if (ost.is_open()) {
290 ost << Tetrahedral_.size() << "\n";
291 ost << "\n";
292
293 for (iter = Tetrahedral_.begin(); iter != Tetrahedral_.end(); ++iter) {
294 Vector3d position;
295 position = (*iter)->getPos();
296
297 ost << "O "
298 << "\t";
299
300 for (unsigned int z = 0; z < position.size(); z++) {
301 ost << position[z] << " "
302 << "\t";
303 }
304 ost << "\n";
305 }
306 ost.close();
307 }
308 }
309 }
310} // namespace OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Real dot(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the dot product of two DynamicVectors.
std::string getPrefix(const std::string &str)