OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
TetrahedralityParamZ.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/TetrahedralityParamZ.hpp"
46
47#include <algorithm>
48#include <fstream>
49#include <vector>
50
51#include "io/DumpReader.hpp"
53#include "utils/simError.h"
54
55using namespace std;
56namespace OpenMD {
57 TetrahedralityParamZ::TetrahedralityParamZ(
58 SimInfo* info, const std::string& filename, const std::string& sele1,
59 const std::string& sele2, double rCut, int nzbins, int axis) :
60 StaticAnalyser(info, filename, nzbins),
61 selectionScript1_(sele1), selectionScript2_(sele2), seleMan1_(info),
62 seleMan2_(info), evaluator1_(info), evaluator2_(info), axis_(axis) {
63 evaluator1_.loadScriptString(sele1);
64 if (!evaluator1_.isDynamic()) {
65 seleMan1_.setSelectionSet(evaluator1_.evaluate());
66 }
67 evaluator2_.loadScriptString(sele2);
68 if (!evaluator2_.isDynamic()) {
69 seleMan2_.setSelectionSet(evaluator2_.evaluate());
70 }
71
72 // Set up cutoff radius:
73 rCut_ = rCut;
74
75 switch (axis_) {
76 case 0:
77 axisLabel_ = "x";
78 break;
79 case 1:
80 axisLabel_ = "y";
81 break;
82 case 2:
83 default:
84 axisLabel_ = "z";
85 break;
86 }
87
88 // fixed number of bins
89 sliceQ_.resize(nBins_);
90 sliceCount_.resize(nBins_);
91 std::fill(sliceQ_.begin(), sliceQ_.end(), 0.0);
92 std::fill(sliceCount_.begin(), sliceCount_.end(), 0);
93
94 setOutputName(getPrefix(filename) + ".Qz");
95 }
96
97 void TetrahedralityParamZ::process() {
98 StuntDouble* sd;
99 StuntDouble* sd2;
100 StuntDouble* sdi;
101 StuntDouble* sdj;
102 int myIndex;
103 Vector3d vec;
104 Vector3d ri, rj, rk, rik, rkj;
105 RealType r;
106 RealType cospsi;
107 RealType Qk;
108 std::vector<std::pair<RealType, StuntDouble*>> myNeighbors;
109 int isd1;
110 int isd2;
111 bool usePeriodicBoundaryConditions_ =
112 info_->getSimParams()->getUsePeriodicBoundaryConditions();
113
114 DumpReader reader(info_, dumpFilename_);
115 int nFrames = reader.getNFrames();
116
117 for (int istep = 0; istep < nFrames; istep += step_) {
118 reader.readFrame(istep);
119 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
120
121 Mat3x3d hmat = currentSnapshot_->getHmat();
122 zBox_.push_back(hmat(axis_, axis_));
123
124 RealType halfBoxZ_ = hmat(axis_, axis_) / 2.0;
125
126 if (evaluator1_.isDynamic()) {
127 seleMan1_.setSelectionSet(evaluator1_.evaluate());
128 }
129
130 if (evaluator2_.isDynamic()) {
131 seleMan2_.setSelectionSet(evaluator2_.evaluate());
132 }
133
134 // outer loop is over the selected StuntDoubles:
135 for (sd = seleMan1_.beginSelected(isd1); sd != NULL;
136 sd = seleMan1_.nextSelected(isd1)) {
137 myIndex = sd->getGlobalIndex();
138
139 Qk = 1.0;
140 myNeighbors.clear();
141
142 for (sd2 = seleMan2_.beginSelected(isd2); sd2 != NULL;
143 sd2 = seleMan2_.nextSelected(isd2)) {
144 if (sd2->getGlobalIndex() != myIndex) {
145 vec = sd->getPos() - sd2->getPos();
146
147 if (usePeriodicBoundaryConditions_)
148 currentSnapshot_->wrapVector(vec);
149
150 r = vec.length();
151
152 // Check to see if neighbor is in bond cutoff
153
154 if (r < rCut_) { myNeighbors.push_back(std::make_pair(r, sd2)); }
155 }
156 }
157
158 // Sort the vector using predicate and std::sort
159 std::sort(myNeighbors.begin(), myNeighbors.end());
160
161 // Use only the 4 closest neighbors to do the rest of the work:
162
163 int nbors = myNeighbors.size() > 4 ? 4 : myNeighbors.size();
164 int nang = int(0.5 * (nbors * (nbors - 1)));
165
166 rk = sd->getPos();
167
168 for (int i = 0; i < nbors - 1; i++) {
169 sdi = myNeighbors[i].second;
170 ri = sdi->getPos();
171 rik = rk - ri;
172 if (usePeriodicBoundaryConditions_) currentSnapshot_->wrapVector(rik);
173
174 rik.normalize();
175
176 for (int j = i + 1; j < nbors; j++) {
177 sdj = myNeighbors[j].second;
178 rj = sdj->getPos();
179 rkj = rk - rj;
180 if (usePeriodicBoundaryConditions_)
181 currentSnapshot_->wrapVector(rkj);
182 rkj.normalize();
183
184 cospsi = dot(rik, rkj);
185
186 // Calculates scaled Qk for each molecule using calculated
187 // angles from 4 or fewer nearest neighbors.
188 Qk -= (pow(cospsi + 1.0 / 3.0, 2) * 2.25 / nang);
189 }
190 }
191
192 if (nang > 0) {
193 if (usePeriodicBoundaryConditions_) currentSnapshot_->wrapVector(rk);
194
195 int binNo =
196 int(nBins_ * (halfBoxZ_ + rk[axis_]) / hmat(axis_, axis_));
197 sliceQ_[binNo] += Qk;
198 sliceCount_[binNo] += 1;
199 }
200 }
201 }
202 writeQz();
203 }
204
205 void TetrahedralityParamZ::writeQz() {
206 // compute average box length:
207
208 RealType zSum = 0.0;
209 for (std::vector<RealType>::iterator j = zBox_.begin(); j != zBox_.end();
210 ++j) {
211 zSum += *j;
212 }
213 RealType zAve = zSum / zBox_.size();
214
215 std::ofstream qZstream(outputFilename_.c_str());
216 if (qZstream.is_open()) {
217 qZstream << "#Tetrahedrality Parameters (" << axisLabel_ << ")\n";
218
219 qZstream << "#nFrames:\t" << zBox_.size() << "\n";
220 qZstream << "#selection 1: (" << selectionScript1_ << ")\n";
221 qZstream << "#selection 2: (" << selectionScript2_ << ")\n";
222 qZstream << "#" << axisLabel_ << "\tQk\n";
223 for (unsigned int i = 0; i < sliceQ_.size(); ++i) {
224 RealType z = zAve * (i + 0.5) / sliceQ_.size();
225 if (sliceCount_[i] != 0) {
226 qZstream << z << "\t" << sliceQ_[i] / sliceCount_[i] << "\n";
227 }
228 }
229
230 } else {
231 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
232 "TetrahedralityParamZ: unable to open %s\n",
233 outputFilename_.c_str());
234 painCave.isFatal = 1;
235 simError();
236 }
237 qZstream.close();
238 }
239} // 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)