OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
GofRAngle2.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/GofRAngle2.hpp"
46
47#include <algorithm>
48#include <fstream>
49#include <sstream>
50
51#include "primitives/Atom.hpp"
52#include "types/MultipoleAdapter.hpp"
53#include "utils/Revision.hpp"
54#include "utils/simError.h"
55
56namespace OpenMD {
57
58 GofRAngle2::GofRAngle2(SimInfo* info, const std::string& filename,
59 const std::string& sele1, const std::string& sele2,
60 RealType len, int nrbins, int nangleBins) :
61 RadialDistrFunc(info, filename, sele1, sele2, nrbins),
62 nAngleBins_(nangleBins), len_(len), doSele3_(false), seleMan3_(info),
63 evaluator3_(info) {
64 setAnalysisType("Radial Distribution Function");
65 setOutputName(getPrefix(filename) + ".grto");
66
67 deltaR_ = len_ / (double)nBins_;
68 deltaCosAngle_ = 2.0 / nAngleBins_;
69
70 std::stringstream params;
71 params << " nBins = " << nBins_ << ", maxLen = " << len_
72 << ", deltaR = " << deltaR_ << ", nAngleBins = " << nAngleBins_
73 << ", deltaCosAngle = " << deltaCosAngle_;
74 const std::string paramString = params.str();
75 setParameterString(paramString);
76
77 histogram_.resize(nBins_);
78 avgGofr_.resize(nBins_);
79 for (unsigned int i = 0; i < nBins_; ++i) {
80 histogram_[i].resize(nAngleBins_);
81 avgGofr_[i].resize(nAngleBins_);
82 for (unsigned int j = 0; j < nAngleBins_; ++j) {
83 histogram_[i][j].resize(nAngleBins_);
84 avgGofr_[i][j].resize(nAngleBins_);
85 }
86 }
87 }
88
89 GofRAngle2::GofRAngle2(SimInfo* info, const std::string& filename,
90 const std::string& sele1, const std::string& sele2,
91 const std::string& sele3, RealType len, int nrbins,
92 int nangleBins) :
93 RadialDistrFunc(info, filename, sele1, sele2, nrbins),
94 nAngleBins_(nangleBins), len_(len), doSele3_(true), seleMan3_(info),
95 evaluator3_(info), selectionScript3_(sele3) {
96 setOutputName(getPrefix(filename) + ".grto");
97
98 deltaR_ = len_ / (double)nBins_;
99 deltaCosAngle_ = 2.0 / nAngleBins_;
100
101 std::stringstream params;
102 params << " nBins = " << nBins_ << ", maxLen = " << len_
103 << ", deltaR = " << deltaR_ << ", nAngleBins = " << nAngleBins_
104 << ", deltaCosAngle = " << deltaCosAngle_;
105 const std::string paramString = params.str();
106 setParameterString(paramString);
107
108 histogram_.resize(nBins_);
109 avgGofr_.resize(nBins_);
110 for (unsigned int i = 0; i < nBins_; ++i) {
111 histogram_[i].resize(nAngleBins_);
112 avgGofr_[i].resize(nAngleBins_);
113 for (unsigned int j = 0; j < nAngleBins_; ++j) {
114 histogram_[i][j].resize(nAngleBins_);
115 avgGofr_[i][j].resize(nAngleBins_);
116 }
117 }
118
119 evaluator3_.loadScriptString(sele3);
120 if (!evaluator3_.isDynamic()) {
121 seleMan3_.setSelectionSet(evaluator3_.evaluate());
122 }
123 }
124
125 void GofRAngle2::processNonOverlapping(SelectionManager& sman1,
126 SelectionManager& sman2) {
127 StuntDouble* sd1;
128 StuntDouble* sd2;
129 StuntDouble* sd3;
130 int i;
131 int j;
132 int k;
133
134 // This is the same as a non-overlapping pairwise loop structure:
135 // for (int i = 0; i < ni ; ++i ) {
136 // for (int j = 0; j < nj; ++j) {}
137 // }
138
139 if (doSele3_) {
140 if (evaluator3_.isDynamic()) {
141 seleMan3_.setSelectionSet(evaluator3_.evaluate());
142 }
143 if (sman1.getSelectionCount() != seleMan3_.getSelectionCount()) {
144 RadialDistrFunc::processNonOverlapping(sman1, sman2);
145 }
146
147 for (sd1 = sman1.beginSelected(i), sd3 = seleMan3_.beginSelected(k);
148 sd1 != NULL && sd3 != NULL;
149 sd1 = sman1.nextSelected(i), sd3 = seleMan3_.nextSelected(k)) {
150 for (sd2 = sman2.beginSelected(j); sd2 != NULL;
151 sd2 = sman2.nextSelected(j)) {
152 collectHistogram(sd1, sd2, sd3);
153 }
154 }
155 } else {
156 RadialDistrFunc::processNonOverlapping(sman1, sman2);
157 }
158 }
159
160 void GofRAngle2::processOverlapping(SelectionManager& sman) {
161 StuntDouble* sd1;
162 StuntDouble* sd2;
163 StuntDouble* sd3;
164 int i;
165 int j;
166 int k;
167
168 // This is the same as a pairwise loop structure:
169 // for (int i = 0; i < n-1 ; ++i ) {
170 // for (int j = i + 1; j < n; ++j) {}
171 // }
172
173 if (doSele3_) {
174 if (evaluator3_.isDynamic()) {
175 seleMan3_.setSelectionSet(evaluator3_.evaluate());
176 }
177 if (sman.getSelectionCount() != seleMan3_.getSelectionCount()) {
178 RadialDistrFunc::processOverlapping(sman);
179 }
180 for (sd1 = sman.beginSelected(i), sd3 = seleMan3_.beginSelected(k);
181 sd1 != NULL && sd3 != NULL;
182 sd1 = sman.nextSelected(i), sd3 = seleMan3_.nextSelected(k)) {
183 for (j = i, sd2 = sman.nextSelected(j); sd2 != NULL;
184 sd2 = sman.nextSelected(j)) {
185 collectHistogram(sd1, sd2, sd3);
186 }
187 }
188 } else {
189 RadialDistrFunc::processOverlapping(sman);
190 }
191 }
192
193 void GofRAngle2::preProcess() {
194 for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
195 for (unsigned int j = 0; j < avgGofr_[i].size(); ++j) {
196 std::fill(avgGofr_[i][j].begin(), avgGofr_[i][j].end(), 0.0);
197 }
198 }
199 }
200
201 void GofRAngle2::initializeHistogram() {
202 for (unsigned int i = 0; i < histogram_.size(); ++i) {
203 for (unsigned int j = 0; j < histogram_.size(); ++j) {
204 std::fill(histogram_[i][j].begin(), histogram_[i][j].end(), 0);
205 }
206 }
207 }
208
209 void GofRAngle2::processHistogram() {
210 int nPairs = getNPairs();
211 RealType volume =
212 info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
213 RealType pairDensity = nPairs / volume;
214 RealType pairConstant = (4.0 * Constants::PI * pairDensity) /
215 (3.0 * (double)nAngleBins_ * (double)nAngleBins_);
216
217 for (unsigned int i = 0; i < histogram_.size(); ++i) {
218 RealType rLower = i * deltaR_;
219 RealType rUpper = rLower + deltaR_;
220 RealType volSlice =
221 (rUpper * rUpper * rUpper) - (rLower * rLower * rLower);
222 RealType nIdeal = volSlice * pairConstant;
223
224 for (unsigned int j = 0; j < histogram_[i].size(); ++j) {
225 for (unsigned int k = 0; k < histogram_[i][j].size(); ++k) {
226 avgGofr_[i][j][k] += RealType(histogram_[i][j][k]) / nIdeal;
227 }
228 }
229 }
230 }
231
232 void GofRAngle2::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
233 bool usePeriodicBoundaryConditions_ =
234 info_->getSimParams()->getUsePeriodicBoundaryConditions();
235
236 if (sd1 == sd2) { return; }
237
238 Vector3d pos1 = sd1->getPos();
239 Vector3d pos2 = sd2->getPos();
240 Vector3d r12 = pos1 - pos2;
241 if (usePeriodicBoundaryConditions_) currentSnapshot_->wrapVector(r12);
242
243 RealType distance = r12.length();
244 int whichRBin = int(distance / deltaR_);
245
246 if (distance <= len_) {
247 if (!sd1->isDirectional()) {
248 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
249 "GofAngle2: attempted to use a non-directional object: %s\n",
250 sd1->getType().c_str());
251 painCave.isFatal = 1;
252 simError();
253 }
254
255 if (!sd2->isDirectional()) {
256 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
257 "GofAngle2: attempted to use a non-directional object: %s\n",
258 sd2->getType().c_str());
259 painCave.isFatal = 1;
260 simError();
261 }
262
263 Vector3d dipole1, dipole2;
264
265 if (sd1->isAtom()) {
266 AtomType* atype1 = static_cast<Atom*>(sd1)->getAtomType();
267 MultipoleAdapter ma1 = MultipoleAdapter(atype1);
268 if (ma1.isDipole())
269 dipole1 = sd1->getDipole();
270 else
271 dipole1 = sd1->getA().transpose() * V3Z;
272
273 } else {
274 dipole1 = sd1->getA().transpose() * V3Z;
275 }
276
277 if (sd2->isAtom()) {
278 AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
279 MultipoleAdapter ma2 = MultipoleAdapter(atype2);
280 if (ma2.isDipole())
281 dipole2 = sd2->getDipole();
282 else
283 dipole2 = sd2->getA().transpose() * V3Z;
284 } else {
285 dipole2 = sd2->getA().transpose() * V3Z;
286 }
287
288 r12.normalize();
289 dipole1.normalize();
290 dipole2.normalize();
291
292 RealType cosAngle1 = dot(r12, dipole1);
293 RealType cosAngle2 = dot(dipole1, dipole2);
294
295 RealType halfBin = (nAngleBins_ - 1) * 0.5;
296 int angleBin1 = int(halfBin * (cosAngle1 + 1.0));
297 int angleBin2 = int(halfBin * (cosAngle2 + 1.0));
298
299 ++histogram_[whichRBin][angleBin1][angleBin2];
300 }
301 }
302
303 void GofRAngle2::collectHistogram(StuntDouble* sd1, StuntDouble* sd2,
304 StuntDouble* sd3) {
305 bool usePeriodicBoundaryConditions_ =
306 info_->getSimParams()->getUsePeriodicBoundaryConditions();
307
308 if (sd1 == sd2) { return; }
309
310 Vector3d p1 = sd1->getPos();
311 Vector3d p3 = sd3->getPos();
312
313 Vector3d c = 0.5 * (p1 + p3);
314 Vector3d r13 = p3 - p1;
315
316 Vector3d r12 = sd2->getPos() - c;
317
318 if (usePeriodicBoundaryConditions_) {
319 currentSnapshot_->wrapVector(r12);
320 currentSnapshot_->wrapVector(r13);
321 }
322
323 RealType distance = r12.length();
324 int whichRBin = int(distance / deltaR_);
325
326 if (distance <= len_) {
327 r12.normalize();
328 r13.normalize();
329
330 if (!sd2->isDirectional()) {
331 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
332 "GofAngle2: attempted to use a non-directional object: %s\n",
333 sd2->getType().c_str());
334 painCave.isFatal = 1;
335 simError();
336 }
337
338 Vector3d dipole2;
339
340 if (sd2->isAtom()) {
341 AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
342 MultipoleAdapter ma2 = MultipoleAdapter(atype2);
343
344 if (ma2.isDipole())
345 dipole2 = sd2->getDipole();
346 else
347 dipole2 = sd2->getA().transpose() * V3Z;
348
349 } else {
350 // Rigid Body:
351 dipole2 = sd2->getA().transpose() * V3Z;
352 }
353
354 dipole2.normalize();
355
356 RealType cosAngle1 = dot(r12, r13);
357 RealType cosAngle2 = dot(r13, dipole2);
358
359 RealType halfBin = (nAngleBins_ - 1) * 0.5;
360 int angleBin1 = int(halfBin * (cosAngle1 + 1.0));
361 int angleBin2 = int(halfBin * (cosAngle2 + 1.0));
362
363 ++histogram_[whichRBin][angleBin1][angleBin2];
364 }
365 }
366
367 void GofRAngle2::writeRdf() {
368 std::ofstream ofs(outputFilename_.c_str());
369 if (ofs.is_open()) {
370 Revision r;
371 ofs << "# " << getAnalysisType() << "\n";
372 ofs << "# OpenMD " << r.getFullRevision() << "\n";
373 ofs << "# " << r.getBuildDate() << "\n";
374 ofs << "# selection script1: \"" << selectionScript1_;
375 ofs << "\"\tselection script2: \"" << selectionScript2_ << "\"";
376 if (doSele3_) {
377 ofs << "\tselection script3: \"" << selectionScript3_ << "\"\n";
378 } else {
379 ofs << "\n";
380 }
381
382 if (!paramString_.empty())
383 ofs << "# parameters: " << paramString_ << "\n";
384
385 for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
386 // RealType r = deltaR_ * (i + 0.5);
387 for (unsigned int j = 0; j < avgGofr_[i].size(); ++j) {
388 // RealType cosAngle1 = -1.0 + (j + 0.5)*deltaCosAngle_;
389 for (unsigned int k = 0; k < avgGofr_[i][j].size(); ++k) {
390 // RealType cosAngle2 = -1.0 + (k + 0.5)*deltaCosAngle_;
391
392 ofs << avgGofr_[i][j][k] / nProcessed_ << "\t";
393 }
394 ofs << "\n";
395 }
396 ofs << "\n";
397 }
398
399 } else {
400 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
401 "GofRAngle2: unable to open %s\n", outputFilename_.c_str());
402 painCave.isFatal = 1;
403 simError();
404 }
405
406 ofs.close();
407 }
408
409} // namespace OpenMD
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)
Real distance(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the distance between two DynamicVectors.