OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
GofRZ.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#include "applications/staticProps/GofRZ.hpp"
46
47#include <algorithm>
48#include <fstream>
49
50#include "utils/simError.h"
51
52namespace OpenMD {
53
54 GofRZ::GofRZ(SimInfo* info, const std::string& filename,
55 const std::string& sele1, const std::string& sele2, RealType len,
56 RealType zlen, int nrbins, int nZBins, int axis) :
57 RadialDistrFunc(info, filename, sele1, sele2, nrbins),
58 len_(len), zLen_(zlen), nZBins_(nZBins), axis_(axis) {
59 setOutputName(getPrefix(filename) + ".gofrz");
60
61 deltaR_ = len_ / (double)nBins_;
62 deltaZ_ = zLen_ / (double)nZBins_;
63
64 histogram_.resize(nBins_);
65 avgGofr_.resize(nBins_);
66 for (unsigned int i = 0; i < nBins_; ++i) {
67 histogram_[i].resize(nZBins_);
68 avgGofr_[i].resize(nZBins_);
69 }
70
71 // Compute complementary axes to the privileged axis
72 xaxis_ = (axis_ + 1) % 3;
73 yaxis_ = (axis_ + 2) % 3;
74
75 // Set the axis label for the privileged axis
76 switch (axis_) {
77 case 0:
78 axisLabel_ = "x";
79 break;
80 case 1:
81 axisLabel_ = "y";
82 break;
83 case 2:
84 default:
85 axisLabel_ = "z";
86 break;
87 }
88 }
89
90 void GofRZ::preProcess() {
91 for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
92 std::fill(avgGofr_[i].begin(), avgGofr_[i].end(), 0);
93 }
94 }
95
96 void GofRZ::initializeHistogram() {
97 npairs_ = 0;
98 for (unsigned int i = 0; i < histogram_.size(); ++i) {
99 std::fill(histogram_[i].begin(), histogram_[i].end(), 0);
100 }
101 }
102
103 void GofRZ::processHistogram() {
104 int nPairs = getNPairs();
105 RealType volume =
106 info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
107 RealType pairDensity = nPairs / volume * 2.0;
108
109 for (unsigned int i = 0; i < histogram_.size(); ++i) {
110 RealType rLower = i * deltaR_;
111 RealType rUpper = rLower + deltaR_;
112 RealType volSlice =
113 Constants::PI * deltaZ_ * ((rUpper * rUpper) - (rLower * rLower));
114 RealType nIdeal = volSlice * pairDensity;
115
116 for (unsigned int j = 0; j < histogram_[i].size(); ++j) {
117 avgGofr_[i][j] += histogram_[i][j] / nIdeal;
118 }
119 }
120 }
121
122 void GofRZ::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
123 if (sd1 == sd2) { return; }
124 bool usePeriodicBoundaryConditions_ =
125 info_->getSimParams()->getUsePeriodicBoundaryConditions();
126
127 Vector3d pos1 = sd1->getPos();
128 Vector3d pos2 = sd2->getPos();
129 Vector3d r12 = pos2 - pos1;
130 if (usePeriodicBoundaryConditions_) currentSnapshot_->wrapVector(r12);
131
132 RealType distance = sqrt(pow(r12[xaxis_], 2) + pow(r12[yaxis_], 2));
133
134 int whichRBin = int(distance / deltaR_);
135
136 if (distance <= len_) {
137 RealType Z = fabs(r12[axis_]);
138
139 if (Z <= zLen_) {
140 int whichZBin = int(Z / deltaZ_);
141
142 ++histogram_[whichRBin][whichZBin];
143 ++npairs_;
144 }
145 }
146 }
147
148 void GofRZ::writeRdf() {
149 std::ofstream rdfStream(outputFilename_.c_str());
150 if (rdfStream.is_open()) {
151 rdfStream << "#radial distribution function\n";
152 rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
153 rdfStream << "selection2: (" << selectionScript2_ << ")\n";
154 rdfStream << "#nBins = " << nBins_ << "\t maxLen = " << len_
155 << "deltaR = " << deltaR_ << "\n";
156 rdfStream << "#n" << axisLabel_ << "Bins =" << nZBins_ << "\t delta"
157 << axisLabel_ << " = " << deltaZ_ << "\n";
158 for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
159 // RealType r = deltaR_ * (i + 0.5);
160
161 for (unsigned int j = 0; j < avgGofr_[i].size(); ++j) {
162 // RealType z = deltaZ_ * (j + 0.5);
163 rdfStream << avgGofr_[i][j] / nProcessed_ << "\t";
164 }
165
166 rdfStream << "\n";
167 }
168
169 } else {
170 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
171 "GofRZ: unable to open %s\n", outputFilename_.c_str());
172 painCave.isFatal = 1;
173 simError();
174 }
175
176 rdfStream.close();
177 }
178
179} // namespace OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
std::string getPrefix(const std::string &str)
Real distance(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the distance between two DynamicVectors.