OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
RCorrFunc.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/dynamicProps/RCorrFunc.hpp"
46
47#include <sstream>
48
49#include "utils/Revision.hpp"
50
51namespace OpenMD {
52 RCorrFunc::RCorrFunc(SimInfo* info, const std::string& filename,
53 const std::string& sele1, const std::string& sele2) :
54 ObjectACF<RealType>(info, filename, sele1, sele2) {
55 setCorrFuncType("Mean Square Displacement");
56 setOutputName(getPrefix(dumpFilename_) + ".rcorr");
57
58 positions_.resize(nFrames_);
59 }
60
61 RCorrFuncZ::RCorrFuncZ(SimInfo* info, const std::string& filename,
62 const std::string& sele1, const std::string& sele2,
63 int nZbins, int axis) :
64 ObjectACF<RealType>(info, filename, sele1, sele2),
65 axis_(axis) {
66 setCorrFuncType("Mean Square Displacement binned by Z");
67 setOutputName(getPrefix(dumpFilename_) + ".rcorrZ");
68
69 positions_.resize(nFrames_);
70 zBins_.resize(nFrames_);
71 nZBins_ = nZbins;
72
73 switch (axis_) {
74 case 0:
75 axisLabel_ = "x";
76 break;
77 case 1:
78 axisLabel_ = "y";
79 break;
80 case 2:
81 default:
82 axisLabel_ = "z";
83 break;
84 }
85
86 std::stringstream params;
87 params << " nzbins = " << nZBins_;
88 const std::string paramString = params.str();
89 setParameterString(paramString);
90
91 histograms_.resize(nTimeBins_);
92 counts_.resize(nTimeBins_);
93
94 idimHistograms_.resize(3);
95 for (unsigned i = 0; i < idimHistograms_.size(); i++) {
96 idimHistograms_[i].resize(nTimeBins_);
97 }
98 for (unsigned int i = 0; i < nTimeBins_; i++) {
99 histograms_[i].resize(nZBins_);
100 counts_[i].resize(nZBins_);
101
102 std::fill(histograms_[i].begin(), histograms_[i].end(), 0.0);
103 std::fill(counts_[i].begin(), counts_[i].end(), 0);
104
105 for (unsigned j = 0; j < 3; j++) {
106 idimHistograms_[j][i].resize(nZBins_);
107 std::fill(idimHistograms_[j][i].begin(), idimHistograms_[j][i].end(),
108 0.0);
109 }
110 }
111 }
112
113 RCorrFuncR::RCorrFuncR(SimInfo* info, const std::string& filename,
114 const std::string& sele1, const std::string& sele2) :
115 ObjectACF<RealType>(info, filename, sele1, sele2) {
116 // Turn on COM calculation in reader:
117 bool ncp = true;
118 reader_->setNeedCOMprops(ncp);
119 setCorrFuncType("MSD (radial projection)");
120 setOutputName(getPrefix(dumpFilename_) + ".r_rcorr");
121 positions_.resize(nFrames_);
122 }
123
124 int RCorrFunc::computeProperty1(int frame, StuntDouble* sd) {
125 positions_[frame].push_back(sd->getPos());
126 return positions_[frame].size() - 1;
127 }
128
129 RealType RCorrFunc::calcCorrVal(int frame1, int frame2, int id1, int id2) {
130 Vector3d diff = positions_[frame2][id2] - positions_[frame1][id1];
131 return diff.lengthSquare();
132 }
133
134 void RCorrFuncZ::computeFrame(int istep) {
135 hmat_ = currentSnapshot_->getHmat();
136 halfBoxZ_ = hmat_(axis_, axis_) / 2.0;
137
138 StuntDouble* sd;
139
140 int isd1, isd2;
141 unsigned int index;
142
143 if (evaluator1_.isDynamic()) {
144 seleMan1_.setSelectionSet(evaluator1_.evaluate());
145 }
146
147 if (uniqueSelections_ && evaluator2_.isDynamic()) {
148 seleMan2_.setSelectionSet(evaluator2_.evaluate());
149 }
150
151 for (sd = seleMan1_.beginSelected(isd1); sd != NULL;
152 sd = seleMan1_.nextSelected(isd1)) {
153 index = computeProperty1(istep, sd);
154 if (index == sele1ToIndex_[istep].size()) {
155 sele1ToIndex_[istep].push_back(sd->getGlobalIndex());
156 } else {
157 sele1ToIndex_[istep].resize(index + 1);
158 sele1ToIndex_[istep][index] = sd->getGlobalIndex();
159 }
160 }
161
162 if (uniqueSelections_) {
163 for (sd = seleMan2_.beginSelected(isd2); sd != NULL;
164 sd = seleMan2_.nextSelected(isd2)) {
165 index = computeProperty1(istep, sd);
166
167 if (index == sele2ToIndex_[istep].size()) {
168 sele2ToIndex_[istep].push_back(sd->getGlobalIndex());
169 } else {
170 sele2ToIndex_[istep].resize(index + 1);
171 sele2ToIndex_[istep][index] = sd->getGlobalIndex();
172 }
173 }
174 }
175 }
176
177 int RCorrFuncZ::computeProperty1(int frame, StuntDouble* sd) {
178 Vector3d pos = sd->getPos();
179 // we need the raw (not wrapped) positions for RMSD:
180 positions_[frame].push_back(sd->getPos());
181
182 if (info_->getSimParams()->getUsePeriodicBoundaryConditions()) {
183 currentSnapshot_->wrapVector(pos);
184 }
185 int zBin = int(nZBins_ * (halfBoxZ_ + pos[axis_]) / hmat_(axis_, axis_));
186 zBins_[frame].push_back(zBin);
187
188 return positions_[frame].size() - 1;
189 }
190
191 void RCorrFuncZ::correlateFrames(int frame1, int frame2, int timeBin) {
192 std::vector<int> s1;
193 std::vector<int> s2;
194
195 std::vector<int>::iterator i1;
196 std::vector<int>::iterator i2;
197
198 s1 = sele1ToIndex_[frame1];
199
200 if (uniqueSelections_)
201 s2 = sele2ToIndex_[frame2];
202 else
203 s2 = sele1ToIndex_[frame2];
204
205 for (i1 = s1.begin(), i2 = s2.begin(); i1 != s1.end() && i2 != s2.end();
206 ++i1, ++i2) {
207 // If the selections are dynamic, they might not have the
208 // same objects in both frames, so we need to roll either of
209 // the selections until we have the same object to
210 // correlate.
211
212 while (i1 != s1.end() && *i1 < *i2) {
213 ++i1;
214 }
215
216 while (i2 != s2.end() && *i2 < *i1) {
217 ++i2;
218 }
219
220 if (i1 == s1.end() || i2 == s2.end()) break;
221
222 calcCorrValImpl(frame1, frame2, i1 - s1.begin(), i2 - s2.begin(),
223 timeBin);
224 }
225 }
226
227 RealType RCorrFuncZ::calcCorrValImpl(int frame1, int frame2, int id1, int id2,
228 int timeBin) {
229 int zBin1 = zBins_[frame1][id1];
230 int zBin2 = zBins_[frame2][id2];
231
232 if (zBin1 == zBin2) {
233 Vector3d diff = positions_[frame2][id2] - positions_[frame1][id1];
234 histograms_[timeBin][zBin1] += diff.lengthSquare();
235
236 for (unsigned i = 0; i < 3; i++) {
237 RealType iDiff =
238 positions_[frame2][id2][i] - positions_[frame1][id1][i];
239 idimHistograms_[i][timeBin][zBin1] += (iDiff * iDiff);
240 }
241
242 counts_[timeBin][zBin1]++;
243 }
244 return 0.0;
245 }
246
247 void RCorrFuncZ::postCorrelate() {
248 for (unsigned int i = 0; i < nTimeBins_; ++i) {
249 for (unsigned int j = 0; j < nZBins_; ++j) {
250 if (counts_[i][j] > 0) {
251 histograms_[i][j] /= counts_[i][j];
252 for (unsigned int k = 0; k < 3; k++) {
253 idimHistograms_[k][i][j] /= counts_[i][j];
254 }
255 } else {
256 histograms_[i][j] = 0;
257 for (unsigned int k = 0; k < 3; k++) {
258 idimHistograms_[k][i][j] = 0;
259 }
260 }
261 }
262 }
263 }
264 void RCorrFuncZ::writeCorrelate() {
265 std::ofstream ofs(getOutputFileName().c_str());
266
267 if (ofs.is_open()) {
268 Revision r;
269
270 ofs << "# " << getCorrFuncType() << "\n";
271 ofs << "# OpenMD " << r.getFullRevision() << "\n";
272 ofs << "# " << r.getBuildDate() << "\n";
273 ofs << "# selection script1: \"" << selectionScript1_;
274 ofs << "\"\tselection script2: \"" << selectionScript2_ << "\"\n";
275 ofs << "# privilegedAxis computed as " << axisLabel_ << " axis \n";
276 if (!paramString_.empty())
277 ofs << "# parameters: " << paramString_ << "\n";
278
279 ofs << "#time\tcorrVal\n";
280
281 for (unsigned int i = 0; i < nTimeBins_; ++i) {
282 ofs << times_[i] - times_[0];
283
284 for (unsigned int j = 0; j < nZBins_; ++j) {
285 ofs << "\t" << histograms_[i][j];
286 }
287 ofs << "\n";
288 }
289
290 ofs << "&\n#time\tcorrValXZ\n";
291
292 for (unsigned int i = 0; i < nTimeBins_; ++i) {
293 ofs << times_[i] - times_[0];
294
295 for (unsigned int j = 0; j < nZBins_; ++j) {
296 ofs << "\t" << idimHistograms_[0][i][j];
297 }
298 ofs << "\n";
299 }
300
301 ofs << "&\n#time\tcorrValYZ\n";
302
303 for (unsigned int i = 0; i < nTimeBins_; ++i) {
304 ofs << times_[i] - times_[0];
305
306 for (unsigned int j = 0; j < nZBins_; ++j) {
307 ofs << "\t" << idimHistograms_[1][i][j];
308 }
309 ofs << "\n";
310 }
311
312 ofs << "&\n#time\tcorrValZZ\n";
313
314 for (unsigned int i = 0; i < nTimeBins_; ++i) {
315 ofs << times_[i] - times_[0];
316
317 for (unsigned int j = 0; j < nZBins_; ++j) {
318 ofs << "\t" << idimHistograms_[2][i][j];
319 }
320 ofs << "\n";
321 }
322
323 } else {
324 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
325 "RCorrFuncZ::writeCorrelate Error: fail to open %s\n",
326 getOutputFileName().c_str());
327 painCave.isFatal = 1;
328 simError();
329 }
330 ofs.close();
331 }
332
333 int RCorrFuncR::computeProperty1(int frame, StuntDouble* sd) {
334 // get the radial vector from the frame's center of mass:
335 Vector3d coord_t = sd->getPos() - sd->getCOM();
336
337 positions_[frame].push_back(coord_t.length());
338 return positions_[frame].size() - 1;
339 }
340
341 RealType RCorrFuncR::calcCorrVal(int frame1, int frame2, int id1, int id2) {
342 RealType dr;
343 dr = positions_[frame2][id2] - positions_[frame1][id1];
344 return dr * dr;
345 }
346} // namespace OpenMD
Real lengthSquare()
Returns the squared length of this vector.
Definition Vector.hpp:399
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
std::string getPrefix(const std::string &str)