OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
SurfaceDiffusion.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/* Surface Diffusion
46 * Attempting to track/measure the surface diffusion rates of particles on...
47 * wait for it.. a surface. This program was initially created to track Platinum
48 * particles moving around a 557 surface. Hence why we are trying to keep the x
49 * and y movement separate.
50 *
51 */
52
53#include "applications/staticProps/SurfaceDiffusion.hpp"
54
55#include <algorithm>
56#include <fstream>
57
58#include "io/DumpReader.hpp"
60#include "utils/simError.h"
61
62namespace OpenMD {
63
64 SurfaceDiffusion::SurfaceDiffusion(SimInfo* info, const std::string& filename,
65 const std::string& sele, RealType) :
66 StaticAnalyser(info, filename, 1),
67 selectionScript_(sele), evaluator_(info), seleMan1_(info) {
68 evaluator_.loadScriptString(sele);
69 if (!evaluator_.isDynamic()) {
70 seleMan1_.setSelectionSet(evaluator_.evaluate());
71 }
72
73 // Depending on the selection 'sele1="select Pt"' need a vector equal to the
74 // number of Platinums in the system (for this specific case)
75 selectionCount_ = seleMan1_.getSelectionCount();
76 cout << "SelectionCount_: " << selectionCount_ << "\n";
77
78 moBool_.resize(selectionCount_);
79 positions_.resize(selectionCount_);
80
81 filename_ = filename;
82 singleMoveDistance_ = 2.0;
83 }
84
85 void SurfaceDiffusion::process() {
86 StuntDouble* sd;
87 bool usePeriodicBoundaryConditions_ =
88 info_->getSimParams()->getUsePeriodicBoundaryConditions();
89
90 DumpReader reader(info_, dumpFilename_);
91 int nFrames = reader.getNFrames();
92 frames_ = 0;
93 nProcessed_ = nFrames / step_;
94
95 // positions_ and moBool_ are 2D arrays, need the second dimension
96 // filled as well
97 for (int i = 0; i < selectionCount_; i++) {
98 moBool_[i].resize(nFrames);
99 positions_[i].resize(nFrames);
100 }
101
102 int iterator;
103 int index = 0;
104 /* Loop over all frames storing the positions in a vec< vec<pos> >
105 * At the end, positions.length() should equal seleMan1_.size() or
106 * w/e And positions[index].length() should equal nFrames (or
107 * nFrames/istep)
108 */
109 for (int istep = 0; istep < nFrames; istep += step_) {
110 frames_++;
111 reader.readFrame(istep);
112 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
113
114 index = 0; // count over atoms since iterators aren't the most
115 // friendly for such plebian things
116 for (sd = seleMan1_.beginSelected(iterator); sd != NULL;
117 sd = seleMan1_.nextSelected(iterator)) {
118 Vector3d pos = sd->getPos();
119 positions_[index][istep] = pos;
120 index++;
121 }
122 }
123
124 cout << "Position Array size: " << positions_.size() << "\n";
125 cout << "Frames analyzed: " << positions_[0].size() << "\n";
126
127 for (std::size_t i = 0; i < positions_.size(); i++) {
128 int frameIndex = positions_[i].size();
129 for (int j = 1; j < frameIndex; j++) {
130 Vector3d posF1 = positions_[i][j - 1];
131 Vector3d posF2 = positions_[i][j];
132 Vector3d diff = posF2 - posF1;
133 if (usePeriodicBoundaryConditions_) {
134 currentSnapshot_->wrapVector(diff);
135 }
136 double dist = diff.length();
137 if (dist > singleMoveDistance_) {
138 moBool_[i][j] = true;
139 } else {
140 moBool_[i][j] = false;
141 }
142 }
143 }
144
145 int mobileAtomCount = 0;
146 for (std::size_t i = 0; i < moBool_.size(); i++) {
147 int frameIndex = moBool_[i].size();
148 bool mobileAtom = false;
149 for (int j = 0; j < frameIndex; j++) {
150 mobileAtom = mobileAtom || moBool_[i][j];
151 }
152 moBool_[i][0] = mobileAtom; // is true if any value later in the
153 // array is true, false otherwise
154 if (mobileAtom) { mobileAtomCount++; }
155 }
156
157 cout << "Mobile atom count: " << mobileAtomCount << "\n";
158
159 // Here I shrink the size of the arrays, why look through 3888,
160 // when you only need ~800. Additionally, all of these are mobile
161 // at some point in time, the others aren't, dead weight and
162 // memory
163 positions2_.resize(mobileAtomCount);
164 moBool2_.resize(mobileAtomCount);
165 int pos2index = 0;
166 for (std::size_t i = 0; i < positions_.size(); i++) {
167 int frameCount = positions_[i].size();
168 if (moBool_[i][0]) {
169 for (int j = 0; j < frameCount; j++) {
170 positions2_[pos2index].push_back(positions_[i][j]);
171 moBool2_[pos2index].push_back(moBool_[i][j]);
172 }
173 pos2index++;
174 }
175 }
176
177 positions_.clear();
178 moBool_.clear();
179
180 cout << "positions_ has been cleared: " << positions_.size() << "\n";
181 cout << "positions2_ has been filled: " << positions2_.size() << "\n";
182 cout << "positions2_ has " << positions2_[0].size() << " frames\n";
183
184 // The important one!
185 positionCorrelation();
186
187 // Write out my data
188 std::ofstream diffStream;
189 setOutputName(getPrefix(filename_) + ".Mdiffusion");
190 diffStream.open(outputFilename_.c_str());
191 diffStream << "#X&Y diffusion amounts\n";
192 diffStream << "#singleMoveDistance_: " << singleMoveDistance_ << "\n";
193 diffStream << "#Number of mobile atoms: " << positions2_.size() << "\n";
194 diffStream << "#time, <x(t)-x(0)>, <y(t)-y(0)>, <r(t)-r(0)>\n";
195
196 for (std::size_t i = 0; i < xHist_.size(); i++) {
197 diffStream << i << ", " << xHist_[i] << ", " << yHist_[i] << ", "
198 << rHist_[i] << "\n";
199 }
200 diffStream.close();
201 }
202
203 void SurfaceDiffusion::positionCorrelation() {
204 RealType xDist = 0.0;
205 RealType yDist = 0.0;
206 RealType rDist = 0.0;
207 int timeShift = 0;
208 Vector3d kPos;
209 Vector3d jPos;
210 // biggest timeShift is positions2_[0].size() - 1?
211 xHist_.clear();
212 yHist_.clear();
213 rHist_.clear();
214 count_.clear();
215 int frameResize = positions2_[0].size();
216 xHist_.resize(frameResize);
217 yHist_.resize(frameResize);
218 rHist_.resize(frameResize);
219 count_.resize(frameResize);
220 // loop over particles
221 // loop over frames starting at j
222 // loop over frames starting at k = j (time shift of 0)
223 for (std::size_t i = 0; i < positions2_.size(); i++) {
224 int frames = positions2_[i].size() - 1; // for counting
225 // properly, otherwise
226 // moBool2_[i][j+1] will
227 // go over
228 for (int j = 0; j < frames; j++) {
229 // if the particle is mobile between j and j + 1, then count
230 // it for all timeShifts
231 if (moBool2_[i][j + 1]) {
232 for (std::size_t k = j; k < positions2_[0].size(); k++) {
233 //<x(t)-x(0)> <y(t)-y(0)> <r(t)-r(0)>
234 // The positions stored are not wrapped, thus I don't need
235 // to worry about pbc
236 // Mean square displacement
237 // So I do want the squared distances
238
239 kPos = positions2_[i][k];
240 jPos = positions2_[i][j];
241 xDist = kPos.x() - jPos.x();
242 xDist = xDist * xDist;
243
244 yDist = kPos.y() - jPos.y();
245 yDist = yDist * yDist;
246
247 rDist = (kPos - jPos).lengthSquare();
248
249 timeShift = k - j;
250 xHist_[timeShift] += xDist;
251 yHist_[timeShift] += yDist;
252 rHist_[timeShift] += rDist;
253 count_[timeShift]++;
254 }
255 }
256 }
257 }
258 cout << "X, Y, R calculated\n";
259
260 for (std::size_t i = 0; i < xHist_.size(); i++) {
261 xHist_[i] = xHist_[i] / (count_[i]);
262 yHist_[i] = yHist_[i] / (count_[i]);
263 rHist_[i] = rHist_[i] / (count_[i]);
264 }
265 cout << "X, Y, R normalized\n";
266 }
267
268} // namespace OpenMD
Real length()
Returns the length of this vector.
Definition Vector.hpp:393
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
std::string getPrefix(const std::string &str)