OpenMD 3.2
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
P2R.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 following paper when you publish your work:
33 *
34 * [1] Drisko et al., J. Open Source Softw. 9, 7004 (2024).
35 *
36 * Good starting points for code and simulation methodology are:
37 *
38 * [2] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
39 * [3] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
40 * [4] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
41 * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
42 * [6] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
43 * [7] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
44 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
45 * [9] Drisko & Gezelter, J. Chem. Theory Comput. 20, 4986-4997 (2024).
46 */
47
48/* Calculates Angle(R) for DirectionalAtoms*/
49
50#include "applications/staticProps/P2R.hpp"
51
52#include <algorithm>
53#include <cmath>
54#include <fstream>
55#include <sstream>
56
57#include "brains/Thermo.hpp"
58#include "io/DumpReader.hpp"
59#include "primitives/Molecule.hpp"
60#include "utils/Revision.hpp"
61#include "utils/simError.h"
62
63namespace OpenMD {
64
65 P2R::P2R(SimInfo* info, const std::string& filename, const std::string& sele1,
66 unsigned int nbins) :
67 StaticAnalyser(info, filename, nbins),
68 doVect_(true), doOffset_(false), selectionScript1_(sele1),
69 seleMan1_(info), seleMan2_(info), evaluator1_(info), evaluator2_(info) {
70 evaluator1_.loadScriptString(sele1);
71 if (!evaluator1_.isDynamic()) {
72 seleMan1_.setSelectionSet(evaluator1_.evaluate());
73 }
74
75 setAnalysisType(
76 "2nd order Legendre Polynomial Correlation using r as reference axis");
77 setOutputName(getPrefix(filename) + ".P2R");
78 std::stringstream params;
79 const std::string paramString = params.str();
80 setParameterString(paramString);
81 }
82
83 P2R::P2R(SimInfo* info, const std::string& filename, const std::string& sele1,
84 const std::string& sele2, unsigned int nbins) :
85 StaticAnalyser(info, filename, nbins),
86 doVect_(false), doOffset_(false), selectionScript1_(sele1),
87 selectionScript2_(sele2), seleMan1_(info), seleMan2_(info),
88 evaluator1_(info), evaluator2_(info) {
89 evaluator1_.loadScriptString(sele1);
90 if (!evaluator1_.isDynamic()) {
91 seleMan1_.setSelectionSet(evaluator1_.evaluate());
92 }
93
94 evaluator2_.loadScriptString(sele2);
95 if (!evaluator2_.isDynamic()) {
96 seleMan2_.setSelectionSet(evaluator2_.evaluate());
97 }
98
99 setAnalysisType(
100 "2nd order Legendre Polynomial Correlation using r as reference axis");
101 setOutputName(getPrefix(filename) + ".P2R");
102 std::stringstream params;
103 const std::string paramString = params.str();
104 setParameterString(paramString);
105 }
106
107 P2R::P2R(SimInfo* info, const std::string& filename, const std::string& sele1,
108 int seleOffset, unsigned int nbins) :
109 StaticAnalyser(info, filename, nbins),
110 doVect_(false), doOffset_(true), doOffset2_(false),
111 selectionScript1_(sele1), seleMan1_(info), seleMan2_(info),
112 evaluator1_(info), evaluator2_(info), seleOffset_(seleOffset) {
113 evaluator1_.loadScriptString(sele1);
114 if (!evaluator1_.isDynamic()) {
115 seleMan1_.setSelectionSet(evaluator1_.evaluate());
116 }
117
118 setAnalysisType(
119 "2nd order Legendre Polynomial Correlation using r as reference axis");
120 setOutputName(getPrefix(filename) + ".P2R");
121 std::stringstream params;
122 const std::string paramString = params.str();
123 setParameterString(paramString);
124 }
125
126 P2R::P2R(SimInfo* info, const std::string& filename, const std::string& sele1,
127 int seleOffset, int seleOffset2, unsigned int nbins) :
128 StaticAnalyser(info, filename, nbins),
129 doVect_(false), doOffset_(true), doOffset2_(true),
130 selectionScript1_(sele1), seleMan1_(info), seleMan2_(info),
131 evaluator1_(info), evaluator2_(info), seleOffset_(seleOffset),
132 seleOffset2_(seleOffset2) {
133 evaluator1_.loadScriptString(sele1);
134 if (!evaluator1_.isDynamic()) {
135 seleMan1_.setSelectionSet(evaluator1_.evaluate());
136 }
137
138 setAnalysisType(
139 "2nd order Legendre Polynomial Correlation using r as reference axis");
140 setOutputName(getPrefix(filename) + ".P2R");
141 std::stringstream params;
142 const std::string paramString = params.str();
143 setParameterString(paramString);
144 }
145
146 void P2R::process() {
147 StuntDouble* sd1;
148 StuntDouble* sd2;
149 int ii;
150 int jj;
151 bool usePeriodicBoundaryConditions_ =
152 info_->getSimParams()->getUsePeriodicBoundaryConditions();
153
154 Thermo thermo(info_);
155 DumpReader reader(info_, dumpFilename_);
156 int nFrames = reader.getNFrames();
157
158 nProcessed_ = nFrames / step_;
159 P2_ = 0.0;
160 count_ = 0;
161
162 for (int istep = 0; istep < nFrames; istep += step_) {
163 reader.readFrame(istep);
164 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
165
166 Vector3d CenterOfMass = thermo.getCom();
167
168 if (evaluator1_.isDynamic()) {
169 seleMan1_.setSelectionSet(evaluator1_.evaluate());
170 }
171
172 if (doVect_) {
173 for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
174 sd1 = seleMan1_.nextSelected(ii)) {
175 Vector3d pos = sd1->getPos();
176
177 Vector3d r1 = CenterOfMass - pos;
178 // only do this if the stunt double actually has a vector associated
179 // with it
180 if (sd1->isDirectional()) {
181 Vector3d vec = sd1->getA().transpose() * V3Z;
182 vec.normalize();
183 r1.normalize();
184 RealType cosangle = dot(r1, vec);
185
186 P2_ += 0.5 * (3.0 * cosangle * cosangle - 1.0);
187 count_++;
188 }
189 }
190 } else {
191 if (doOffset_) {
192 for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
193 sd1 = seleMan1_.nextSelected(ii)) {
194 // This will require careful rewriting if StaticProps is
195 // ever parallelized. For an example, see
196 // Thermo::getTaggedAtomPairDistance
197 Vector3d r1;
198
199 if (doOffset2_) {
200 int sd1Aind = sd1->getGlobalIndex() + seleOffset2_;
201 StuntDouble* sd1A = info_->getIOIndexToIntegrableObject(sd1Aind);
202 r1 = CenterOfMass - sd1A->getPos();
203 } else {
204 r1 = CenterOfMass - sd1->getPos();
205 }
206
207 if (usePeriodicBoundaryConditions_)
208 currentSnapshot_->wrapVector(r1);
209
210 int sd2Index = sd1->getGlobalIndex() + seleOffset_;
211 sd2 = info_->getIOIndexToIntegrableObject(sd2Index);
212
213 Vector3d r2 = CenterOfMass - sd2->getPos();
214 if (usePeriodicBoundaryConditions_)
215 currentSnapshot_->wrapVector(r2);
216
217 Vector3d rc = 0.5 * (r1 + r2);
218 if (usePeriodicBoundaryConditions_)
219 currentSnapshot_->wrapVector(rc);
220
221 Vector3d vec = r1 - r2;
222 if (usePeriodicBoundaryConditions_)
223 currentSnapshot_->wrapVector(vec);
224
225 rc.normalize();
226 vec.normalize();
227 RealType cosangle = dot(rc, vec);
228 P2_ += 0.5 * (3.0 * cosangle * cosangle - 1.0);
229 count_++;
230 }
231 } else {
232 if (evaluator2_.isDynamic()) {
233 seleMan2_.setSelectionSet(evaluator2_.evaluate());
234 }
235
236 if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount()) {
237 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
238 "In frame %d, the number of selected StuntDoubles are\n"
239 "\tnot the same in --sele1 and sele2\n",
240 istep);
241 painCave.severity = OPENMD_INFO;
242 painCave.isFatal = 0;
243 simError();
244 }
245
246 for (sd1 = seleMan1_.beginSelected(ii),
247 sd2 = seleMan2_.beginSelected(jj);
248 sd1 != NULL && sd2 != NULL; sd1 = seleMan1_.nextSelected(ii),
249 sd2 = seleMan2_.nextSelected(jj)) {
250 Vector3d r1 = CenterOfMass - sd1->getPos();
251 if (usePeriodicBoundaryConditions_)
252 currentSnapshot_->wrapVector(r1);
253
254 Vector3d r2 = CenterOfMass - sd2->getPos();
255 if (usePeriodicBoundaryConditions_)
256 currentSnapshot_->wrapVector(r2);
257
258 Vector3d rc = 0.5 * (r1 + r2);
259 if (usePeriodicBoundaryConditions_)
260 currentSnapshot_->wrapVector(rc);
261
262 Vector3d vec = r1 - r2;
263 if (usePeriodicBoundaryConditions_)
264 currentSnapshot_->wrapVector(vec);
265
266 rc.normalize();
267 vec.normalize();
268 RealType cosangle = dot(rc, vec);
269 P2_ += 0.5 * (3.0 * cosangle * cosangle - 1.0);
270 count_++;
271 }
272 }
273 }
274 }
275 processHistogram();
276 writeP2R();
277 }
278
279 void P2R::processHistogram() {
280 if (count_ > 0) P2_ /= count_;
281 }
282
283 void P2R::writeP2R() {
284 std::ofstream ofs(outputFilename_.c_str());
285 if (ofs.is_open()) {
286 Revision rev;
287
288 ofs << "# " << getAnalysisType() << "\n";
289 ofs << "# OpenMD " << rev.getFullRevision() << "\n";
290 ofs << "# " << rev.getBuildDate() << "\n";
291 ofs << "#nFrames:\t" << nProcessed_ << "\n";
292 ofs << "#selection1: (" << selectionScript1_ << ")";
293 if (!doVect_) { ofs << "\tselection2: (" << selectionScript2_ << ")"; }
294 ofs << "\n";
295 if (!paramString_.empty())
296 ofs << "# parameters: " << paramString_ << "\n";
297
298 ofs << "#P2 correlation:\n";
299 ofs << P2_ << "\n";
300
301 } else {
302 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
303 "P2R: unable to open %s\n", outputFilename_.c_str());
304 painCave.isFatal = 1;
305 simError();
306 }
307
308 ofs.close();
309 }
310
311 P2Z::P2Z(SimInfo* info, const std::string& filename, const std::string& sele1,
312 unsigned int nbins, int axis) :
313 P2R(info, filename, sele1, nbins),
314 axis_(axis) {
315 switch (axis_) {
316 case 0:
317 axisLabel_ = "x";
318 break;
319 case 1:
320 axisLabel_ = "y";
321 break;
322 case 2:
323 default:
324 axisLabel_ = "z";
325 break;
326 }
327
328 setAnalysisType("2nd order Legendre Polynomial Correlation using " +
329 axisLabel_ + " as reference axis");
330 setOutputName(getPrefix(filename) + ".P2Z");
331 std::stringstream params;
332 const std::string paramString = params.str();
333 setParameterString(paramString);
334 }
335
336 P2Z::P2Z(SimInfo* info, const std::string& filename, const std::string& sele1,
337 const std::string& sele2, unsigned int nbins, int axis) :
338 P2R(info, filename, sele1, sele2, nbins),
339 axis_(axis) {
340 switch (axis_) {
341 case 0:
342 axisLabel_ = "x";
343 break;
344 case 1:
345 axisLabel_ = "y";
346 break;
347 case 2:
348 default:
349 axisLabel_ = "z";
350 break;
351 }
352
353 setAnalysisType("2nd order Legendre Polynomial Correlation using " +
354 axisLabel_ + " as reference axis");
355 setOutputName(getPrefix(filename) + ".P2Z");
356 std::stringstream params;
357 const std::string paramString = params.str();
358 setParameterString(paramString);
359 }
360
361 P2Z::P2Z(SimInfo* info, const std::string& filename, const std::string& sele1,
362 int seleOffset, unsigned int nbins, int axis) :
363 P2R(info, filename, sele1, seleOffset, nbins),
364 axis_(axis) {
365 switch (axis_) {
366 case 0:
367 axisLabel_ = "x";
368 break;
369 case 1:
370 axisLabel_ = "y";
371 break;
372 case 2:
373 default:
374 axisLabel_ = "z";
375 break;
376 }
377
378 setAnalysisType("2nd order Legendre Polynomial Correlation using " +
379 axisLabel_ + " as reference axis");
380 setOutputName(getPrefix(filename) + ".P2Z");
381 std::stringstream params;
382 const std::string paramString = params.str();
383 setParameterString(paramString);
384 }
385
386 P2Z::P2Z(SimInfo* info, const std::string& filename, const std::string& sele1,
387 int seleOffset, int seleOffset2, unsigned int nbins, int axis) :
388 P2R(info, filename, sele1, seleOffset, seleOffset2, nbins),
389 axis_(axis) {
390 switch (axis_) {
391 case 0:
392 axisLabel_ = "x";
393 break;
394 case 1:
395 axisLabel_ = "y";
396 break;
397 case 2:
398 default:
399 axisLabel_ = "z";
400 break;
401 }
402
403 setAnalysisType("2nd order Legendre Polynomial Correlation using " +
404 axisLabel_ + " as reference axis");
405 setOutputName(getPrefix(filename) + ".P2Z");
406 std::stringstream params;
407 const std::string paramString = params.str();
408 setParameterString(paramString);
409 }
410
411 void P2Z::process() {
412 Vector3d zhat = V3Zero;
413 zhat[axis_] = 1.0; // unit vector along preferred axis
414
415 StuntDouble* sd1;
416 StuntDouble* sd2;
417 int ii;
418 int jj;
419 bool usePeriodicBoundaryConditions_ =
420 info_->getSimParams()->getUsePeriodicBoundaryConditions();
421
422 DumpReader reader(info_, dumpFilename_);
423 int nFrames = reader.getNFrames();
424
425 nProcessed_ = nFrames / step_;
426 P2_ = 0.0;
427 count_ = 0;
428
429 for (int istep = 0; istep < nFrames; istep += step_) {
430 reader.readFrame(istep);
431 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
432
433 if (evaluator1_.isDynamic()) {
434 seleMan1_.setSelectionSet(evaluator1_.evaluate());
435 }
436
437 if (doVect_) {
438 for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
439 sd1 = seleMan1_.nextSelected(ii)) {
440 // only do this if the stunt double actually has a vector associated
441 // with it
442 if (sd1->isDirectional()) {
443 Vector3d vec = sd1->getA().transpose() * V3Z;
444 vec.normalize();
445 RealType cosangle = dot(zhat, vec);
446
447 P2_ += 0.5 * (3.0 * cosangle * cosangle - 1.0);
448 count_++;
449 }
450 }
451 } else {
452 if (doOffset_) {
453 for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
454 sd1 = seleMan1_.nextSelected(ii)) {
455 // This will require careful rewriting if StaticProps is
456 // ever parallelized. For an example, see
457 // Thermo::getTaggedAtomPairDistance
458 Vector3d r1;
459
460 if (doOffset2_) {
461 int sd1Aind = sd1->getGlobalIndex() + seleOffset2_;
462 StuntDouble* sd1A = info_->getIOIndexToIntegrableObject(sd1Aind);
463 r1 = sd1A->getPos();
464 } else {
465 r1 = sd1->getPos();
466 }
467
468 int sd2Index = sd1->getGlobalIndex() + seleOffset_;
469 sd2 = info_->getIOIndexToIntegrableObject(sd2Index);
470
471 Vector3d r2 = sd2->getPos();
472 Vector3d vec = r1 - r2;
473
474 if (usePeriodicBoundaryConditions_)
475 currentSnapshot_->wrapVector(vec);
476
477 vec.normalize();
478 RealType cosangle = dot(zhat, vec);
479 P2_ += 0.5 * (3.0 * cosangle * cosangle - 1.0);
480 count_++;
481 }
482 } else {
483 if (evaluator2_.isDynamic()) {
484 seleMan2_.setSelectionSet(evaluator2_.evaluate());
485 }
486
487 if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount()) {
488 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
489 "In frame %d, the number of selected StuntDoubles are\n"
490 "\tnot the same in --sele1 and sele2\n",
491 istep);
492 painCave.severity = OPENMD_INFO;
493 painCave.isFatal = 0;
494 simError();
495 }
496
497 for (sd1 = seleMan1_.beginSelected(ii),
498 sd2 = seleMan2_.beginSelected(jj);
499 sd1 != NULL && sd2 != NULL; sd1 = seleMan1_.nextSelected(ii),
500 sd2 = seleMan2_.nextSelected(jj)) {
501 Vector3d r1 = sd1->getPos();
502 Vector3d r2 = sd2->getPos();
503 Vector3d vec = r1 - r2;
504
505 if (usePeriodicBoundaryConditions_)
506 currentSnapshot_->wrapVector(vec);
507
508 vec.normalize();
509 RealType cosangle = dot(zhat, vec);
510 P2_ += 0.5 * (3.0 * cosangle * cosangle - 1.0);
511 count_++;
512 }
513 }
514 }
515 }
516 processHistogram();
517 writeP2R();
518 }
519
520} // namespace OpenMD
OpenMD::SimInfo
One of the heavy-weight classes of OpenMD, SimInfo maintains objects and variables relating to the cu...
Definition SimInfo.hpp:96
OpenMD::Vector::normalize
void normalize()
Normalizes this vector in place.
Definition Vector.hpp:406
OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Definition ActionCorrFunc.cpp:63
OpenMD::dot
Real dot(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the dot product of two DynamicVectors.
Definition DynamicVector.hpp:452
OpenMD::getPrefix
std::string getPrefix(const std::string &str)
Definition StringUtils.cpp:207