OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
PotDiff.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
46
47#include <algorithm>
48#include <functional>
49
52#include "io/DumpReader.hpp"
54#include "types/FixedChargeAdapter.hpp"
55#include "types/FluctuatingChargeAdapter.hpp"
56#include "utils/simError.h"
57
58namespace OpenMD {
59
60 PotDiff::PotDiff(SimInfo* info, const std::string& filename,
61 const std::string& sele) :
62 StaticAnalyser(info, filename, 1),
63 selectionScript_(sele), seleMan_(info), evaluator_(info) {
64 StuntDouble* sd;
65 int i;
66
67 setOutputName(getPrefix(filename) + ".potDiff");
68
69 // The PotDiff is computed by negating the charge on the atom type
70 // using fluctuating charge values. If we don't have any
71 // fluctuating charges in the simulation, we need to expand
72 // storage to hold them.
73 int atomStorageLayout = info_->getAtomStorageLayout();
74 int rigidBodyStorageLayout = info->getRigidBodyStorageLayout();
75 int cutoffGroupStorageLayout = info->getCutoffGroupStorageLayout();
76
77 atomStorageLayout |= DataStorage::dslFlucQPosition;
78 atomStorageLayout |= DataStorage::dslFlucQVelocity;
79 atomStorageLayout |= DataStorage::dslFlucQForce;
80
81 info_->setAtomStorageLayout(atomStorageLayout);
82 info_->setSnapshotManager(new SimSnapshotManager(info_, atomStorageLayout,
83 rigidBodyStorageLayout,
84 cutoffGroupStorageLayout));
85
86 // now we have to figure out which AtomTypes to convert to fluctuating
87 // charges
88 evaluator_.loadScriptString(sele);
89 seleMan_.setSelectionSet(evaluator_.evaluate());
90 for (sd = seleMan_.beginSelected(i); sd != NULL;
91 sd = seleMan_.nextSelected(i)) {
92 AtomType* at = static_cast<Atom*>(sd)->getAtomType();
94 if (fqa.isFluctuatingCharge()) {
95 selectionWasFlucQ_.push_back(true);
96 } else {
97 selectionWasFlucQ_.push_back(false);
98 // make a fictitious fluctuating charge with an unphysical
99 // charge mass and slaterN, but we need to zero out the
100 // electronegativity and hardness to remove the self
101 // contribution:
102 fqa.makeFluctuatingCharge(1.0e9, 0.0, 0.0, 1);
103 sd->setFlucQPos(0.0);
104 }
105 }
106 info_->getSnapshotManager()->advance();
107 }
108
110 StuntDouble* sd;
111 int j;
112
113 diff_.clear();
114 DumpReader reader(info_, dumpFilename_);
115 int nFrames = reader.getNFrames();
116
117 // We'll need the force manager to compute the potential
118 ForceManager* forceMan = new ForceManager(info_);
119
120 // We'll need thermo to report the potential
121 Thermo* thermo = new Thermo(info_);
122
123 for (int i = 0; i < nFrames; i += step_) {
124 reader.readFrame(i);
125 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
126
127 for (sd = seleMan_.beginSelected(j); sd != NULL;
128 sd = seleMan_.nextSelected(j)) {
129 if (!selectionWasFlucQ_[j]) { sd->setFlucQPos(0.0); }
130 }
131
132 forceMan->calcForces();
133 RealType pot1 = thermo->getPotential();
134
135 if (evaluator_.isDynamic()) {
136 seleMan_.setSelectionSet(evaluator_.evaluate());
137 }
138
139 for (sd = seleMan_.beginSelected(j); sd != NULL;
140 sd = seleMan_.nextSelected(j)) {
141 AtomType* at = static_cast<Atom*>(sd)->getAtomType();
142
145
146 RealType charge = 0.0;
147
148 if (fca.isFixedCharge()) charge += fca.getCharge();
149 if (fqa.isFluctuatingCharge()) charge += sd->getFlucQPos();
150
151 sd->setFlucQPos(-charge);
152 }
153
154 currentSnapshot_->clearDerivedProperties();
155 forceMan->calcForces();
156 RealType pot2 = thermo->getPotential();
157 RealType diff = pot2 - pot1;
158
159 data_.add(diff);
160 diff_.push_back(diff);
161 times_.push_back(currentSnapshot_->getTime());
162
163 info_->getSnapshotManager()->advance();
164 }
165
166 writeDiff();
167 }
168
169 void PotDiff::writeDiff() {
170 std::ofstream ofs(outputFilename_.c_str(), std::ios::binary);
171 if (ofs.is_open()) {
172 RealType mu = data_.getAverage();
173 RealType sigma = data_.getStdDev();
174 RealType m95 = data_.get95percentConfidenceInterval();
175
176 ofs << "#potDiff\n";
177 ofs << "#selection: (" << selectionScript_ << ")\n";
178 ofs << "# <diff> = " << mu << "\n";
179 ofs << "# StdDev = " << sigma << "\n";
180 ofs << "# 95% confidence interval = " << m95 << "\n";
181 ofs << "# t\tdiff[t]\n";
182 for (unsigned int i = 0; i < diff_.size(); ++i) {
183 ofs << times_[i] << "\t" << diff_[i] << "\n";
184 }
185
186 } else {
187 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
188 "PotDiff: unable to open %s\n", outputFilename_.c_str());
189 painCave.isFatal = 1;
190 simError();
191 }
192 ofs.close();
193 }
194
195} // namespace OpenMD
StaticAnalyser for Potential Energy changes with charges turned off.
void getAverage(ResultType &ret)
compute the Mean
void get95percentConfidenceInterval(ResultType &ret)
return the 95% confidence interval:
virtual void add(ElementType const &val)
Accumulate another value.
void getStdDev(ResultType &ret)
compute error of average value
AtomType is what OpenMD looks to for unchanging data about an atom.
Definition AtomType.hpp:66
ForceManager is responsible for calculating both the short range (bonded) interactions and long range...
virtual void process()
Process the data.
Definition PotDiff.cpp:109
PotDiff(SimInfo *info, const std::string &filename, const std::string &sele)
Default constructor.
Definition PotDiff.cpp:60
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
void setAtomStorageLayout(int asl)
Sets the storage layouts (computed by SimCreator)
Definition SimInfo.hpp:256
int getAtomStorageLayout()
Returns the storage layouts (computed by SimCreator)
Definition SimInfo.hpp:251
void setSnapshotManager(SnapshotManager *sman)
Sets the snapshot manager.
Definition SimInfo.cpp:971
SnapshotManager * getSnapshotManager()
Returns the snapshot manager.
Definition SimInfo.hpp:248
"brains/SimSnapshotManager.hpp"
void clearDerivedProperties()
sets the state of the computed properties to false
Definition Snapshot.cpp:135
Snapshot * getCurrentSnapshot()
Returns the pointer of current snapshot.
"Don't move, or you're dead! Stand up! Captain, we've got them!"
void setFlucQPos(RealType charge)
Sets the current fluctuating charge of this stuntDouble.
RealType getFlucQPos()
Returns the current fluctuating charge of this stuntDouble.
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
std::string getPrefix(const std::string &str)