OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
LegendreCorrFuncZ.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/dynamicProps/LegendreCorrFuncZ.hpp"
46
47#include <sstream>
48
50#include "utils/Revision.hpp"
51#include "utils/simError.h"
52
53namespace OpenMD {
54 LegendreCorrFuncZ::LegendreCorrFuncZ(SimInfo* info,
55 const std::string& filename,
56 const std::string& sele1,
57 const std::string& sele2, int order,
58 int nZbins, int axis) :
59 ObjectACF<Vector3d>(info, filename, sele1, sele2),
60 nZBins_(nZbins), axis_(axis) {
61 setCorrFuncType("Legendre Correlation Function of Z");
62 setOutputName(getPrefix(dumpFilename_) + ".lcorrZ");
63
64 std::stringstream params;
65 params << " order = " << order << ", nzbins = " << nZBins_;
66 const std::string paramString = params.str();
67 setParameterString(paramString);
68
69 if (!uniqueSelections_) { seleMan2_ = seleMan1_; }
70
71 // Compute complementary axes to the privileged axis
72 xaxis_ = (axis_ + 1) % 3;
73 yaxis_ = (axis_ + 2) % 3;
74
75 switch (axis_) {
76 case 0:
77 axisLabel_ = "x";
78 break;
79 case 1:
80 axisLabel_ = "y";
81 break;
82 case 2:
83 default:
84 axisLabel_ = "z";
85 break;
86 }
87
88 rotMats_.resize(nTimeBins_);
89 zbin_.resize(nTimeBins_);
90 histogram_.resize(nTimeBins_);
91 counts_.resize(nTimeBins_);
92 for (unsigned int i = 0; i < nTimeBins_; i++) {
93 histogram_[i].resize(nZBins_);
94 std::fill(histogram_[i].begin(), histogram_[i].end(), 0.0);
95 counts_[i].resize(nZBins_);
96 std::fill(counts_[i].begin(), counts_[i].end(), 0);
97 }
98 LegendrePolynomial polynomial(order);
99 legendre_ = polynomial.getLegendrePolynomial(order);
100 }
101
102 void LegendreCorrFuncZ::computeFrame(int frame) {
103 Mat3x3d hmat = currentSnapshot_->getHmat();
104 boxZ_ = hmat(axis_, axis_);
105 halfBoxZ_ = boxZ_ / 2.0;
106
108 }
109
110 int LegendreCorrFuncZ::computeProperty1(int frame, StuntDouble* sd) {
111 RotMat3x3d A = sd->getA();
112 rotMats_[frame].push_back(A);
113
114 Vector3d pos = sd->getPos();
115 if (info_->getSimParams()->getUsePeriodicBoundaryConditions())
116 currentSnapshot_->wrapVector(pos);
117 int zBin = int(nZBins_ * (halfBoxZ_ + pos[axis_]) / boxZ_);
118 zbin_[frame].push_back(zBin);
119
120 return rotMats_[frame].size() - 1;
121 }
122
123 Vector3d LegendreCorrFuncZ::calcCorrVal(int frame1, int frame2, int id1,
124 int id2) {
125 Vector3d v1x = rotMats_[frame1][id1].getRow(xaxis_);
126 Vector3d v1y = rotMats_[frame1][id1].getRow(yaxis_);
127 Vector3d v1z = rotMats_[frame1][id1].getRow(axis_);
128
129 Vector3d v2x = rotMats_[frame2][id2].getRow(xaxis_);
130 Vector3d v2y = rotMats_[frame2][id2].getRow(yaxis_);
131 Vector3d v2z = rotMats_[frame2][id2].getRow(axis_);
132
133 RealType uxprod =
134 legendre_.evaluate(dot(v1x, v2x) / (v1x.length() * v2x.length()));
135 RealType uyprod =
136 legendre_.evaluate(dot(v1y, v2y) / (v1y.length() * v2y.length()));
137 RealType uzprod =
138 legendre_.evaluate(dot(v1z, v2z) / (v1z.length() * v2z.length()));
139
140 return Vector3d(uxprod, uyprod, uzprod);
141 }
142
143 void LegendreCorrFuncZ::correlateFrames(int frame1, int frame2, int timeBin) {
144 std::vector<int> s1;
145 std::vector<int> s2;
146
147 std::vector<int>::iterator i1;
148 std::vector<int>::iterator i2;
149
150 Vector3d corrVal(0.0);
151
152 s1 = sele1ToIndex_[frame1];
153
154 if (uniqueSelections_)
155 s2 = sele2ToIndex_[frame2];
156 else
157 s2 = sele1ToIndex_[frame2];
158
159 for (i1 = s1.begin(), i2 = s2.begin(); i1 != s1.end() && i2 != s2.end();
160 ++i1, ++i2) {
161 // If the selections are dynamic, they might not have the
162 // same objects in both frames, so we need to roll either of
163 // the selections until we have the same object to
164 // correlate.
165
166 while (i1 != s1.end() && *i1 < *i2) {
167 ++i1;
168 }
169
170 while (i2 != s2.end() && *i2 < *i1) {
171 ++i2;
172 }
173
174 if (i1 == s1.end() || i2 == s2.end()) break;
175
176 corrVal = calcCorrVal(frame1, frame2, i1 - s1.begin(), i2 - s2.begin());
177 int zBin = zbin_[frame1][i1 - s1.begin()];
178 histogram_[timeBin][zBin] += corrVal;
179 counts_[timeBin][zBin]++;
180 }
181 }
182
183 void LegendreCorrFuncZ::postCorrelate() {
184 for (unsigned int i = 0; i < nTimeBins_; ++i) {
185 for (unsigned int j = 0; j < nZBins_; ++j) {
186 if (counts_[i][j] > 0) { histogram_[i][j] /= counts_[i][j]; }
187 }
188 }
189 }
190
191 void LegendreCorrFuncZ::validateSelection(SelectionManager&) {
192 StuntDouble* sd;
193 int i;
194 for (sd = seleMan1_.beginSelected(i); sd != NULL;
195 sd = seleMan1_.nextSelected(i)) {
196 if (!sd->isDirectionalAtom()) {
197 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
198 "LegendreCorrFunc::validateSelection Error: "
199 "selected atoms are not Directional\n");
200 painCave.isFatal = 1;
201 simError();
202 }
203 }
204 }
205
206 void LegendreCorrFuncZ::writeCorrelate() {
207 std::ofstream ofs(getOutputFileName().c_str());
208
209 if (ofs.is_open()) {
210 Revision r;
211
212 ofs << "# " << getCorrFuncType() << "\n";
213 ofs << "# OpenMD " << r.getFullRevision() << "\n";
214 ofs << "# " << r.getBuildDate() << "\n";
215 ofs << "# selection script1: \"" << selectionScript1_;
216 ofs << "\"\tselection script2: \"" << selectionScript2_ << "\"\n";
217 ofs << "# privilegedAxis computed as " << axisLabel_ << " axis \n";
218 if (!paramString_.empty())
219 ofs << "# parameters: " << paramString_ << "\n";
220
221 ofs << "#time\tPn(costheta_z)\n";
222
223 for (unsigned int i = 0; i < nTimeBins_; ++i) {
224 ofs << times_[i] - times_[0];
225
226 for (unsigned int j = 0; j < nZBins_; ++j) {
227 ofs << "\t" << histogram_[i][j](2);
228 }
229 ofs << "\n";
230 }
231
232 } else {
233 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
234 "LegendreCorrFuncZ::writeCorrelate Error: failed to open %s\n",
235 getOutputFileName().c_str());
236 painCave.isFatal = 1;
237 simError();
238 }
239 ofs.close();
240 }
241} // namespace OpenMD
A collection of Legendre Polynomials.
Real evaluate(const Real &x)
Calculates the value of this Polynomial evaluated at the given x value.
StuntDouble * nextSelected(int &i)
Finds the next selected StuntDouble in the selection.
StuntDouble * beginSelected(int &i)
Finds the first selected StuntDouble in the selection.
One of the heavy-weight classes of OpenMD, SimInfo maintains objects and variables relating to the cu...
Definition SimInfo.hpp:93
Mat3x3d getHmat()
Returns the H-Matrix.
Definition Snapshot.cpp:214
void wrapVector(Vector3d &v)
Wrapping the vector according to periodic boundary condition.
Definition Snapshot.cpp:337
"Don't move, or you're dead! Stand up! Captain, we've got them!"
RotMat3x3d getA()
Returns the current rotation matrix of this stuntDouble.
bool isDirectionalAtom()
Tests if this stuntDouble is an directional atom.
Vector3d getPos()
Returns the current position of this stuntDouble.
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)