OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
CurrentDensity.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/*
46 * Computes the current density for the selected atom
47 * Created by Cody R. Drisko on 06/14/19.
48 */
49
50#include "applications/staticProps/CurrentDensity.hpp"
51
52#include <algorithm>
53#include <fstream>
54#include <string>
55#include <vector>
56
57#include "brains/Thermo.hpp"
58#include "io/DumpReader.hpp"
60#include "types/FixedChargeAdapter.hpp"
61#include "types/FluctuatingChargeAdapter.hpp"
62#include "utils/StringUtils.hpp"
63#include "utils/simError.h"
64
65namespace OpenMD {
66
67 CurrentDensity::CurrentDensity(SimInfo* info, const std::string& filename,
68 const std::string& sele, int nbins, int axis) :
69 StaticAnalyser(info, filename, nbins),
70 selectionScript_(sele), evaluator_(info), seleMan_(info), thermo_(info),
71 axis_(axis) {
72 evaluator_.loadScriptString(sele);
73 if (!evaluator_.isDynamic()) {
74 seleMan_.setSelectionSet(evaluator_.evaluate());
75 }
76
77 // fixed number of bins
78 sliceSDLists_.resize(nBins_);
79 currentDensity_.resize(nBins_);
80
81 switch (axis_) {
82 case 0:
83 axisLabel_ = "x";
84 break;
85 case 1:
86 axisLabel_ = "y";
87 break;
88 case 2:
89 default:
90 axisLabel_ = "z";
91 break;
92 }
93
94 setOutputName(getPrefix(filename) + ".Jc");
95 }
96
97 void CurrentDensity::process() {
98 StuntDouble* sd;
99 int ii;
100
101 bool usePeriodicBoundaryConditions_ =
102 info_->getSimParams()->getUsePeriodicBoundaryConditions();
103
104 DumpReader reader(info_, dumpFilename_);
105 int nFrames = reader.getNFrames();
106 nProcessed_ = nFrames / step_;
107 overallCurrentDensity_ = 0;
108
109 for (int istep = 0; istep < nFrames; istep += step_) {
110 reader.readFrame(istep);
111 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
112 Vector3d COMvel = thermo_.getComVel();
113
114 for (unsigned int i = 0; i < nBins_; i++) {
115 sliceSDLists_[i].clear();
116 }
117
118 RealType sliceVolume = currentSnapshot_->getVolume() / nBins_;
119 Mat3x3d hmat = currentSnapshot_->getHmat();
120 zBox_.push_back(hmat(axis_, axis_));
121
122 if (evaluator_.isDynamic()) {
123 seleMan_.setSelectionSet(evaluator_.evaluate());
124 }
125
126 // determine which atom belongs to which slice
127 for (sd = seleMan_.beginSelected(ii); sd != NULL;
128 sd = seleMan_.nextSelected(ii)) {
129 int binNo;
130 Vector3d pos = sd->getPos();
131
132 if (usePeriodicBoundaryConditions_) {
133 currentSnapshot_->wrapVector(pos);
134 binNo =
135 int(nBins_ * (pos[axis_] / hmat(axis_, axis_) + 0.5)) % nBins_;
136 sliceSDLists_[binNo].push_back(sd);
137 }
138 }
139
140 // loop over the slices to calculate the densities
141 for (unsigned int i = 0; i < nBins_; i++) {
142 RealType binJc = 0;
143
144 for (unsigned int j = 0; j < sliceSDLists_[i].size(); ++j) {
145 RealType q = 0.0;
146 Atom* atom = static_cast<Atom*>(sliceSDLists_[i][j]);
147
148 AtomType* atomType = atom->getAtomType();
149
150 if (sliceSDLists_[i][j]->isAtom()) {
152 if (fca.isFixedCharge()) q = fca.getCharge();
153
155 if (fqa.isFluctuatingCharge()) q += atom->getFlucQPos();
156
157 Vector3d vel = sliceSDLists_[i][j]->getVel();
158 binJc += q * (vel[axis_] - COMvel[axis_]);
159 }
160 }
161
162 // Units of (e / Ang^2 / fs)
163 currentDensity_[i] += binJc / sliceVolume;
164 overallCurrentDensity_ += currentDensity_[i];
165 }
166 }
167
168 writeCurrentDensity();
169 }
170
171 void CurrentDensity::writeCurrentDensity() {
172 // compute average box length:
173 std::vector<RealType>::iterator j;
174 RealType zSum = 0.0;
175 for (j = zBox_.begin(); j != zBox_.end(); ++j) {
176 zSum += *j;
177 }
178 RealType zAve = zSum / zBox_.size();
179
180 std::ofstream rdfStream(outputFilename_.c_str());
181 if (rdfStream.is_open()) {
182 rdfStream << "#Current Density = "
183 << overallCurrentDensity_ / (nBins_ * nProcessed_)
184 << " e / Ang^2 / fs.\n";
185 rdfStream << "#J_c(" << axisLabel_ << ")\n";
186 rdfStream << "#nFrames:\t" << nProcessed_ << "\n";
187 rdfStream << "#selection: (" << selectionScript_ << ")\n";
188 rdfStream << "#" << axisLabel_ << "\tcurrent density\n";
189
190 for (unsigned int i = 0; i < currentDensity_.size(); ++i) {
191 RealType z = zAve * (i + 0.5) / currentDensity_.size();
192 rdfStream << z << "\t" << currentDensity_[i] / nProcessed_ << "\n";
193 }
194
195 } else {
196 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
197 "CurrentDensity: unable to open %s\n", outputFilename_.c_str());
198 painCave.isFatal = 1;
199 simError();
200 }
201
202 rdfStream.close();
203 }
204} // namespace OpenMD
AtomType is what OpenMD looks to for unchanging data about an atom.
Definition AtomType.hpp:66
bool isDynamic()
Tests if the result from evaluation of script is dynamic.
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
SnapshotManager * getSnapshotManager()
Returns the snapshot manager.
Definition SimInfo.hpp:248
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
Snapshot * getCurrentSnapshot()
Returns the pointer of current snapshot.
"Don't move, or you're dead! Stand up! Captain, we've got them!"
Vector3d getPos()
Returns the current position of this stuntDouble.
Vector3d getComVel()
Returns the velocity of center of mass of the whole system.
Definition Thermo.cpp:775
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
std::string getPrefix(const std::string &str)