OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
pAngle.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 Rho(theta) */
46
47#include "applications/staticProps/pAngle.hpp"
48
49#include <algorithm>
50#include <fstream>
51
52#include "brains/Thermo.hpp"
53#include "io/DumpReader.hpp"
55#include "utils/simError.h"
56
57namespace OpenMD {
58
59 pAngle::pAngle(SimInfo* info, const std::string& filename,
60 const std::string& sele1, int nthetabins) :
61 StaticAnalyser(info, filename, nthetabins),
62 doVect_(true), doOffset_(false), selectionScript1_(sele1),
63 seleMan1_(info), seleMan2_(info), evaluator1_(info), evaluator2_(info),
64 nThetaBins_(nthetabins) {
65 setOutputName(getPrefix(filename) + ".pAngle");
66
67 evaluator1_.loadScriptString(sele1);
68 if (!evaluator1_.isDynamic()) {
69 seleMan1_.setSelectionSet(evaluator1_.evaluate());
70 }
71
72 count_.resize(nThetaBins_);
73 histogram_.resize(nThetaBins_);
74 }
75
76 pAngle::pAngle(SimInfo* info, const std::string& filename,
77 const std::string& sele1, const std::string& sele2,
78 int nthetabins) :
79 StaticAnalyser(info, filename, nthetabins),
80 doVect_(false), doOffset_(false), selectionScript1_(sele1),
81 selectionScript2_(sele2), seleMan1_(info), seleMan2_(info),
82 evaluator1_(info), evaluator2_(info), nThetaBins_(nthetabins) {
83 setOutputName(getPrefix(filename) + ".pAngle");
84
85 evaluator1_.loadScriptString(sele1);
86 if (!evaluator1_.isDynamic()) {
87 seleMan1_.setSelectionSet(evaluator1_.evaluate());
88 }
89
90 evaluator2_.loadScriptString(sele2);
91 if (!evaluator2_.isDynamic()) {
92 seleMan2_.setSelectionSet(evaluator2_.evaluate());
93 }
94
95 count_.resize(nThetaBins_);
96 histogram_.resize(nThetaBins_);
97 }
98
99 pAngle::pAngle(SimInfo* info, const std::string& filename,
100 const std::string& sele1, int seleOffset, int nthetabins) :
101 StaticAnalyser(info, filename, nthetabins),
102 doVect_(false), doOffset_(true), doOffset2_(false),
103 selectionScript1_(sele1), seleMan1_(info), seleMan2_(info),
104 evaluator1_(info), evaluator2_(info), seleOffset_(seleOffset),
105 nThetaBins_(nthetabins) {
106 setOutputName(getPrefix(filename) + ".pAngle");
107
108 evaluator1_.loadScriptString(sele1);
109 if (!evaluator1_.isDynamic()) {
110 seleMan1_.setSelectionSet(evaluator1_.evaluate());
111 }
112
113 count_.resize(nThetaBins_);
114 histogram_.resize(nThetaBins_);
115 }
116
117 pAngle::pAngle(SimInfo* info, const std::string& filename,
118 const std::string& sele1, int seleOffset, int seleOffset2,
119 int nthetabins) :
120 StaticAnalyser(info, filename, nthetabins),
121 doVect_(false), doOffset_(true), doOffset2_(true),
122 selectionScript1_(sele1), seleMan1_(info), seleMan2_(info),
123 evaluator1_(info), evaluator2_(info), seleOffset_(seleOffset),
124 seleOffset2_(seleOffset2), nThetaBins_(nthetabins) {
125 setOutputName(getPrefix(filename) + ".pAngle");
126
127 evaluator1_.loadScriptString(sele1);
128 if (!evaluator1_.isDynamic()) {
129 seleMan1_.setSelectionSet(evaluator1_.evaluate());
130 }
131
132 count_.resize(nThetaBins_);
133 histogram_.resize(nThetaBins_);
134 }
135
136 void pAngle::process() {
137 StuntDouble* sd1;
138 StuntDouble* sd2;
139 int ii;
140 int jj;
141 bool usePeriodicBoundaryConditions_ =
142 info_->getSimParams()->getUsePeriodicBoundaryConditions();
143
144 Thermo thermo(info_);
145 DumpReader reader(info_, dumpFilename_);
146 int nFrames = reader.getNFrames();
147
148 nProcessed_ = nFrames / step_;
149
150 std::fill(histogram_.begin(), histogram_.end(), 0.0);
151 std::fill(count_.begin(), count_.end(), 0);
152
153 for (int istep = 0; istep < nFrames; istep += step_) {
154 reader.readFrame(istep);
155 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
156
157 Vector3d CenterOfMass = thermo.getCom();
158
159 if (evaluator1_.isDynamic()) {
160 seleMan1_.setSelectionSet(evaluator1_.evaluate());
161 }
162
163 if (doVect_) {
164 for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
165 sd1 = seleMan1_.nextSelected(ii)) {
166 Vector3d pos = sd1->getPos();
167
168 Vector3d r1 = CenterOfMass - pos;
169 // only do this if the stunt double actually has a vector associated
170 // with it
171 if (sd1->isDirectional()) {
172 Vector3d vec = sd1->getA().transpose() * V3Z;
173 vec.normalize();
174 r1.normalize();
175 RealType cosangle = dot(r1, vec);
176
177 int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
178 count_[binNo]++;
179 }
180 }
181 } else {
182 if (doOffset_) {
183 for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
184 sd1 = seleMan1_.nextSelected(ii)) {
185 // This will require careful rewriting if StaticProps is
186 // ever parallelized. For an example, see
187 // Thermo::getTaggedAtomPairDistance
188 Vector3d r1;
189
190 if (doOffset2_) {
191 int sd1Aind = sd1->getGlobalIndex() + seleOffset2_;
192 StuntDouble* sd1A = info_->getIOIndexToIntegrableObject(sd1Aind);
193 r1 = CenterOfMass - sd1A->getPos();
194 } else {
195 r1 = CenterOfMass - sd1->getPos();
196 }
197
198 if (usePeriodicBoundaryConditions_)
199 currentSnapshot_->wrapVector(r1);
200
201 int sd2Index = sd1->getGlobalIndex() + seleOffset_;
202 sd2 = info_->getIOIndexToIntegrableObject(sd2Index);
203
204 Vector3d r2 = CenterOfMass - sd2->getPos();
205 if (usePeriodicBoundaryConditions_)
206 currentSnapshot_->wrapVector(r1);
207
208 Vector3d rc = 0.5 * (r1 + r2);
209 if (usePeriodicBoundaryConditions_)
210 currentSnapshot_->wrapVector(rc);
211
212 Vector3d vec = r1 - r2;
213 if (usePeriodicBoundaryConditions_)
214 currentSnapshot_->wrapVector(vec);
215
216 rc.normalize();
217 vec.normalize();
218 RealType cosangle = dot(rc, vec);
219 int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
220 count_[binNo]++;
221 }
222 } else {
223 if (evaluator2_.isDynamic()) {
224 seleMan2_.setSelectionSet(evaluator2_.evaluate());
225 }
226
227 if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount()) {
228 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
229 "In frame %d, the number of selected StuntDoubles are\n"
230 "\tnot the same in --sele1 and sele2\n",
231 istep);
232 painCave.severity = OPENMD_INFO;
233 painCave.isFatal = 0;
234 simError();
235 }
236
237 for (sd1 = seleMan1_.beginSelected(ii),
238 sd2 = seleMan2_.beginSelected(jj);
239 sd1 != NULL && sd2 != NULL; sd1 = seleMan1_.nextSelected(ii),
240 sd2 = seleMan2_.nextSelected(jj)) {
241 Vector3d r1 = CenterOfMass - sd1->getPos();
242 if (usePeriodicBoundaryConditions_)
243 currentSnapshot_->wrapVector(r1);
244
245 Vector3d r2 = CenterOfMass - sd2->getPos();
246 if (usePeriodicBoundaryConditions_)
247 currentSnapshot_->wrapVector(r1);
248
249 Vector3d rc = 0.5 * (r1 + r2);
250 if (usePeriodicBoundaryConditions_)
251 currentSnapshot_->wrapVector(rc);
252
253 Vector3d vec = r1 - r2;
254 if (usePeriodicBoundaryConditions_)
255 currentSnapshot_->wrapVector(vec);
256
257 rc.normalize();
258 vec.normalize();
259 RealType cosangle = dot(rc, vec);
260 int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
261 count_[binNo]++;
262 }
263 }
264 }
265 }
266
267 processHistogram();
268 writeProbs();
269 }
270
271 void pAngle::processHistogram() {
272 int atot = 0;
273 for (unsigned int i = 0; i < count_.size(); ++i)
274 atot += count_[i];
275
276 for (unsigned int i = 0; i < count_.size(); ++i) {
277 histogram_[i] = double(count_[i] / double(atot));
278 }
279 }
280
281 void pAngle::writeProbs() {
282 std::ofstream rdfStream(outputFilename_.c_str());
283 if (rdfStream.is_open()) {
284 rdfStream << "#pAngle\n";
285 rdfStream << "#nFrames:\t" << nProcessed_ << "\n";
286 rdfStream << "#selection1: (" << selectionScript1_ << ")";
287 if (!doVect_) {
288 rdfStream << "\tselection2: (" << selectionScript2_ << ")";
289 }
290 rdfStream << "\n";
291 rdfStream << "#cos(theta)\tp(cos(theta))\n";
292 RealType dct = 2.0 / histogram_.size();
293 for (unsigned int i = 0; i < histogram_.size(); ++i) {
294 RealType ct = -1.0 + (2.0 * i + 1) / (histogram_.size());
295 rdfStream << ct << "\t" << histogram_[i] / dct << "\n";
296 }
297
298 } else {
299 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
300 "pAngle: unable to open %s\n", outputFilename_.c_str());
301 painCave.isFatal = 1;
302 simError();
303 }
304
305 rdfStream.close();
306 }
307
308} // 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)