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