OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
PipeDensity.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 PipeDensity, rho(axis2,axis3) in the box */
46
47#include "applications/staticProps/PipeDensity.hpp"
48
49#include <algorithm>
50#include <fstream>
51
52#include "io/DumpReader.hpp"
54#include "utils/simError.h"
55
56namespace OpenMD {
57
58 PipeDensity::PipeDensity(SimInfo* info, const std::string& filename,
59 const std::string& sele, int nbins, int nbins2,
60 int axis) :
61 StaticAnalyser(info, filename, nbins2),
62 selectionScript_(sele), evaluator_(info), seleMan_(info), nBins2_(nbins),
63 axis_(axis) {
64 evaluator_.loadScriptString(sele);
65 if (!evaluator_.isDynamic()) {
66 seleMan_.setSelectionSet(evaluator_.evaluate());
67 }
68
69 // fixed number of bins
70
71 sliceSDLists_.resize(nBins2_);
72 density_.resize(nBins2_);
73 for (unsigned int i = 0; i < nBins2_; ++i) {
74 sliceSDLists_[i].resize(nBins_);
75 density_[i].resize(nBins_);
76 }
77
78 // Compute complementary axes to the privileged axis
79 axis1_ = (axis_ + 1) % 3;
80 axis2_ = (axis_ + 2) % 3;
81
82 // Set the axis labels for the non-privileged axes
83 switch (axis_) {
84 case 0:
85 axisLabel1_ = "y";
86 axisLabel2_ = "z";
87 break;
88 case 1:
89 axisLabel1_ = "z";
90 axisLabel2_ = "x";
91 break;
92 case 2:
93 default:
94 axisLabel1_ = "x";
95 axisLabel2_ = "y";
96 break;
97 }
98
99 setOutputName(getPrefix(filename) + ".PipeDensity");
100 }
101
102 void PipeDensity::process() {
103 StuntDouble* sd;
104 int ii;
105
106 bool usePeriodicBoundaryConditions_ =
107 info_->getSimParams()->getUsePeriodicBoundaryConditions();
108
109 DumpReader reader(info_, dumpFilename_);
110 int nFrames = reader.getNFrames();
111 nProcessed_ = nFrames / step_;
112
113 for (int istep = 0; istep < nFrames; istep += step_) {
114 reader.readFrame(istep);
115 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
116
117 for (unsigned int i = 0; i < nBins2_; i++) {
118 for (unsigned int j = 0; j < nBins_; j++) {
119 sliceSDLists_[i][j].clear();
120 }
121 }
122
123 RealType sliceVolume = currentSnapshot_->getVolume() / (nBins2_ * nBins_);
124 Mat3x3d hmat = currentSnapshot_->getHmat();
125
126 RealType halfBox1_ = hmat(axis1_, axis1_) / 2.0;
127 RealType halfBox2_ = hmat(axis2_, axis2_) / 2.0;
128
129 if (evaluator_.isDynamic()) {
130 seleMan_.setSelectionSet(evaluator_.evaluate());
131 }
132
133 // wrap the stuntdoubles into a cell
134 for (sd = seleMan_.beginSelected(ii); sd != NULL;
135 sd = seleMan_.nextSelected(ii)) {
136 Vector3d pos = sd->getPos();
137 if (usePeriodicBoundaryConditions_) currentSnapshot_->wrapVector(pos);
138 sd->setPos(pos);
139 }
140
141 // determine which atom belongs to which slice
142 for (sd = seleMan_.beginSelected(ii); sd != NULL;
143 sd = seleMan_.nextSelected(ii)) {
144 Vector3d pos = sd->getPos();
145 // shift molecules by half a box to have bins start at 0
146 int binNo1 =
147 int(nBins2_ * (halfBox1_ + pos[axis1_]) / hmat(axis1_, axis1_));
148 int binNo2 =
149 int(nBins_ * (halfBox2_ + pos[axis2_]) / hmat(axis2_, axis2_));
150 sliceSDLists_[binNo1][binNo2].push_back(sd);
151 }
152
153 // loop over the slices to calculate the densities
154 for (unsigned int i = 0; i < nBins2_; i++) {
155 for (unsigned int j = 0; j < nBins_; j++) {
156 RealType totalMass = 0;
157 for (unsigned int k = 0; k < sliceSDLists_[i][j].size(); ++k) {
158 totalMass += sliceSDLists_[i][j][k]->getMass();
159 }
160 density_[i][j] += totalMass / sliceVolume;
161 }
162 }
163 }
164
165 writeDensity();
166 }
167
168 void PipeDensity::writeDensity() {
169 std::ofstream rdfStream(outputFilename_.c_str());
170
171 if (rdfStream.is_open()) {
172 rdfStream << "#PipeDensity\n";
173 rdfStream << "#nFrames:\t" << nProcessed_ << "\n";
174 rdfStream << "#selection: (" << selectionScript_ << ")\n";
175 rdfStream << "#density (" << axisLabel1_ << "," << axisLabel2_ << ")\n";
176 for (unsigned int i = 0; i < density_.size(); ++i) {
177 for (unsigned int j = 0; j < density_[i].size(); ++j) {
178 rdfStream << Constants::densityConvert * density_[i][j] / nProcessed_;
179 rdfStream << "\t";
180 }
181 rdfStream << "\n";
182 }
183
184 } else {
185 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
186 "PipeDensity: unable to open %s\n", outputFilename_.c_str());
187 painCave.isFatal = 1;
188 simError();
189 }
190
191 rdfStream.close();
192 }
193} // namespace OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
std::string getPrefix(const std::string &str)