OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
AngleR.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/* Calculates Angle(R) for DirectionalAtoms*/
46
47#include "applications/staticProps/AngleR.hpp"
48
49#include <algorithm>
50#include <cmath>
51#include <fstream>
52#include <sstream>
53
54#include "brains/Thermo.hpp"
55#include "io/DumpReader.hpp"
57#include "utils/Revision.hpp"
58#include "utils/simError.h"
59
60namespace OpenMD {
61
62 AngleR::AngleR(SimInfo* info, const std::string& filename,
63 const std::string& sele1, RealType len, int nrbins) :
64 StaticAnalyser(info, filename, nrbins),
65 doVect_(true), doOffset_(false), selectionScript1_(sele1),
66 seleMan1_(info), seleMan2_(info), evaluator1_(info), evaluator2_(info),
67 len_(len), nRBins_(nrbins) {
68 evaluator1_.loadScriptString(sele1);
69 if (!evaluator1_.isDynamic()) {
70 seleMan1_.setSelectionSet(evaluator1_.evaluate());
71 }
72
73 deltaR_ = len_ / nRBins_;
74
75 histogram_.resize(nRBins_);
76 count_.resize(nRBins_);
77 avgAngleR_.resize(nRBins_);
78
79 setAnalysisType("radial density function Angle(r)");
80 setOutputName(getPrefix(filename) + ".AngleR");
81 std::stringstream params;
82 params << " len = " << len_ << ", nrbins = " << nRBins_;
83 const std::string paramString = params.str();
84 setParameterString(paramString);
85 }
86
87 AngleR::AngleR(SimInfo* info, const std::string& filename,
88 const std::string& sele1, const std::string& sele2,
89 RealType len, int nrbins) :
90 StaticAnalyser(info, filename, nrbins),
91 doVect_(false), doOffset_(false), selectionScript1_(sele1),
92 selectionScript2_(sele2), seleMan1_(info), seleMan2_(info),
93 evaluator1_(info), evaluator2_(info), len_(len), nRBins_(nrbins) {
94 setOutputName(getPrefix(filename) + ".AngleR");
95
96 evaluator1_.loadScriptString(sele1);
97 if (!evaluator1_.isDynamic()) {
98 seleMan1_.setSelectionSet(evaluator1_.evaluate());
99 }
100
101 evaluator2_.loadScriptString(sele2);
102 if (!evaluator2_.isDynamic()) {
103 seleMan2_.setSelectionSet(evaluator2_.evaluate());
104 }
105
106 deltaR_ = len_ / nRBins_;
107
108 histogram_.resize(nRBins_);
109 count_.resize(nRBins_);
110 avgAngleR_.resize(nRBins_);
111
112 setAnalysisType("radial density function Angle(r)");
113 setOutputName(getPrefix(filename) + ".AngleR");
114 std::stringstream params;
115 params << " len = " << len_ << ", nrbins = " << nRBins_;
116 const std::string paramString = params.str();
117 setParameterString(paramString);
118 }
119
120 AngleR::AngleR(SimInfo* info, const std::string& filename,
121 const std::string& sele1, int seleOffset, RealType len,
122 int nrbins) :
123 StaticAnalyser(info, filename, nrbins),
124 doVect_(false), doOffset_(true), doOffset2_(false),
125 selectionScript1_(sele1), seleMan1_(info), seleMan2_(info),
126 evaluator1_(info), evaluator2_(info), seleOffset_(seleOffset), len_(len),
127 nRBins_(nrbins) {
128 setOutputName(getPrefix(filename) + ".AngleR");
129
130 evaluator1_.loadScriptString(sele1);
131 if (!evaluator1_.isDynamic()) {
132 seleMan1_.setSelectionSet(evaluator1_.evaluate());
133 }
134
135 deltaR_ = len_ / nRBins_;
136
137 histogram_.resize(nRBins_);
138 count_.resize(nRBins_);
139 avgAngleR_.resize(nRBins_);
140
141 setAnalysisType("radial density function Angle(r)");
142 setOutputName(getPrefix(filename) + ".AngleR");
143 std::stringstream params;
144 params << " len = " << len_ << ", nrbins = " << nRBins_;
145 const std::string paramString = params.str();
146 setParameterString(paramString);
147 }
148
149 AngleR::AngleR(SimInfo* info, const std::string& filename,
150 const std::string& sele1, int seleOffset, int seleOffset2,
151 RealType len, int nrbins) :
152 StaticAnalyser(info, filename, nrbins),
153 doVect_(false), doOffset_(true), doOffset2_(true),
154 selectionScript1_(sele1), seleMan1_(info), seleMan2_(info),
155 evaluator1_(info), evaluator2_(info), seleOffset_(seleOffset),
156 seleOffset2_(seleOffset2), len_(len), nRBins_(nrbins) {
157 setOutputName(getPrefix(filename) + ".AngleR");
158
159 evaluator1_.loadScriptString(sele1);
160 if (!evaluator1_.isDynamic()) {
161 seleMan1_.setSelectionSet(evaluator1_.evaluate());
162 }
163
164 histogram_.resize(nRBins_);
165 count_.resize(nRBins_);
166 avgAngleR_.resize(nRBins_);
167
168 setAnalysisType("radial density function Angle(r)");
169 setOutputName(getPrefix(filename) + ".AngleR");
170 std::stringstream params;
171 params << " len = " << len_ << ", nrbins = " << nRBins_;
172 const std::string paramString = params.str();
173 setParameterString(paramString);
174 }
175
176 void AngleR::process() {
177 StuntDouble* sd1;
178 StuntDouble* sd2;
179 int ii;
180 int jj;
181 RealType distance;
182 bool usePeriodicBoundaryConditions_ =
183 info_->getSimParams()->getUsePeriodicBoundaryConditions();
184
185 Thermo thermo(info_);
186 DumpReader reader(info_, dumpFilename_);
187 int nFrames = reader.getNFrames();
188
189 nProcessed_ = nFrames / step_;
190
191 std::fill(avgAngleR_.begin(), avgAngleR_.end(), 0.0);
192 std::fill(histogram_.begin(), histogram_.end(), 0.0);
193 std::fill(count_.begin(), count_.end(), 0);
194
195 for (int istep = 0; istep < nFrames; istep += step_) {
196 reader.readFrame(istep);
197 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
198
199 Vector3d CenterOfMass = thermo.getCom();
200
201 if (evaluator1_.isDynamic()) {
202 seleMan1_.setSelectionSet(evaluator1_.evaluate());
203 }
204
205 if (doVect_) {
206 for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
207 sd1 = seleMan1_.nextSelected(ii)) {
208 Vector3d pos = sd1->getPos();
209
210 Vector3d r1 = CenterOfMass - pos;
211 distance = r1.length();
212 // only do this if the stunt double actually has a vector associated
213 // with it
214 if (sd1->isDirectional()) {
215 Vector3d vec = sd1->getA().transpose() * V3Z;
216 vec.normalize();
217 r1.normalize();
218 RealType cosangle = dot(r1, vec);
219
220 if (distance < len_) {
221 int whichBin = int(distance / deltaR_);
222 histogram_[whichBin] += cosangle;
223 count_[whichBin] += 1;
224 }
225 }
226 }
227 } else {
228 if (doOffset_) {
229 for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
230 sd1 = seleMan1_.nextSelected(ii)) {
231 // This will require careful rewriting if StaticProps is
232 // ever parallelized. For an example, see
233 // Thermo::getTaggedAtomPairDistance
234 Vector3d r1;
235
236 if (doOffset2_) {
237 int sd1Aind = sd1->getGlobalIndex() + seleOffset2_;
238 StuntDouble* sd1A = info_->getIOIndexToIntegrableObject(sd1Aind);
239 r1 = CenterOfMass - sd1A->getPos();
240 } else {
241 r1 = CenterOfMass - sd1->getPos();
242 }
243
244 if (usePeriodicBoundaryConditions_)
245 currentSnapshot_->wrapVector(r1);
246
247 int sd2Index = sd1->getGlobalIndex() + seleOffset_;
248 sd2 = info_->getIOIndexToIntegrableObject(sd2Index);
249
250 Vector3d r2 = CenterOfMass - sd2->getPos();
251 if (usePeriodicBoundaryConditions_)
252 currentSnapshot_->wrapVector(r1);
253
254 Vector3d rc = 0.5 * (r1 + r2);
255 if (usePeriodicBoundaryConditions_)
256 currentSnapshot_->wrapVector(rc);
257
258 distance = rc.length();
259
260 Vector3d vec = r1 - r2;
261 if (usePeriodicBoundaryConditions_)
262 currentSnapshot_->wrapVector(vec);
263
264 rc.normalize();
265 vec.normalize();
266 RealType cosangle = dot(rc, vec);
267 if (distance < len_) {
268 int whichBin = int(distance / deltaR_);
269 histogram_[whichBin] += cosangle;
270 count_[whichBin] += 1;
271 }
272 }
273 } else {
274 if (evaluator2_.isDynamic()) {
275 seleMan2_.setSelectionSet(evaluator2_.evaluate());
276 }
277
278 if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount()) {
279 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
280 "In frame %d, the number of selected StuntDoubles are\n"
281 "\tnot the same in --sele1 and sele2\n",
282 istep);
283 painCave.severity = OPENMD_INFO;
284 painCave.isFatal = 0;
285 simError();
286 }
287
288 for (sd1 = seleMan1_.beginSelected(ii),
289 sd2 = seleMan2_.beginSelected(jj);
290 sd1 != NULL && sd2 != NULL; sd1 = seleMan1_.nextSelected(ii),
291 sd2 = seleMan2_.nextSelected(jj)) {
292 Vector3d r1 = CenterOfMass - sd1->getPos();
293 if (usePeriodicBoundaryConditions_)
294 currentSnapshot_->wrapVector(r1);
295
296 Vector3d r2 = CenterOfMass - sd2->getPos();
297 if (usePeriodicBoundaryConditions_)
298 currentSnapshot_->wrapVector(r1);
299
300 Vector3d rc = 0.5 * (r1 + r2);
301 if (usePeriodicBoundaryConditions_)
302 currentSnapshot_->wrapVector(rc);
303
304 Vector3d vec = r1 - r2;
305 if (usePeriodicBoundaryConditions_)
306 currentSnapshot_->wrapVector(vec);
307
308 distance = rc.length();
309
310 rc.normalize();
311 vec.normalize();
312 RealType cosangle = dot(rc, vec);
313 if (distance < len_) {
314 int whichBin = int(distance / deltaR_);
315 histogram_[whichBin] += cosangle;
316 count_[whichBin] += 1;
317 }
318 }
319 }
320 }
321 }
322 processHistogram();
323 writeAngleR();
324 }
325
326 void AngleR::processHistogram() {
327 for (unsigned int i = 0; i < histogram_.size(); ++i) {
328 if (count_[i] > 0)
329 avgAngleR_[i] += histogram_[i] / count_[i];
330 else
331 avgAngleR_[i] = 0.0;
332
333 std::cerr << " count = " << count_[i] << " avgAngle = " << avgAngleR_[i]
334 << "\n";
335 }
336 }
337
338 void AngleR::writeAngleR() {
339 std::ofstream ofs(outputFilename_.c_str());
340 if (ofs.is_open()) {
341 Revision rev;
342
343 ofs << "# " << getAnalysisType() << "\n";
344 ofs << "# OpenMD " << rev.getFullRevision() << "\n";
345 ofs << "# " << rev.getBuildDate() << "\n";
346 ofs << "#nFrames:\t" << nProcessed_ << "\n";
347 ofs << "#selection1: (" << selectionScript1_ << ")";
348 if (!doVect_) { ofs << "\tselection2: (" << selectionScript2_ << ")"; }
349 if (!paramString_.empty())
350 ofs << "# parameters: " << paramString_ << "\n";
351
352 ofs << "#r\tcorrValue\n";
353 for (unsigned int i = 0; i < avgAngleR_.size(); ++i) {
354 RealType r = deltaR_ * (i + 0.5);
355 ofs << r << "\t" << avgAngleR_[i] << "\n";
356 }
357
358 } else {
359 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
360 "AngleR: unable to open %s\n", outputFilename_.c_str());
361 painCave.isFatal = 1;
362 simError();
363 }
364
365 ofs.close();
366 }
367
368} // namespace OpenMD
void normalize()
Normalizes this vector in place.
Definition Vector.hpp:402
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)
Real distance(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the distance between two DynamicVectors.