OpenMD 3.1
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
NitrileFrequencyMap.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/NitrileFrequencyMap.hpp"
46
47#include <algorithm>
48#include <fstream>
49
50#include "brains/Thermo.hpp"
51#include "io/DumpReader.hpp"
53#include "utils/simError.h"
54
55namespace OpenMD {
56
57 NitrileFrequencyMap::NitrileFrequencyMap(SimInfo* info,
58 const std::string& filename,
59 const std::string& sele1,
60 int nbins) :
61 StaticAnalyser(info, filename, nbins),
62 info_(info), selectionScript1_(sele1), seleMan1_(info_),
63 evaluator1_(info_) {
64 setOutputName(getPrefix(filename) + ".freqs");
65
66 evaluator1_.loadScriptString(sele1);
67 if (!evaluator1_.isDynamic()) {
68 seleMan1_.setSelectionSet(evaluator1_.evaluate());
69 }
70
71 count_.resize(nBins_);
72 histogram_.resize(nBins_);
73
74 freqs_.resize(info_->getNGlobalMolecules());
75
76 minFreq_ = -50;
77 maxFreq_ = 50;
78
79 // Values from Choi et. al. "Nitrile and thiocyanate IR probes:
80 // Quantum chemistry calculation studies and multivariate
81 // least-square ļ¬tting analysis," J. Chem. Phys. 128, 134506 (2008).
82 //
83 // These map site electrostatic potentials onto frequency shifts
84 // in the same energy units that one computes the total potential.
85
86 frequencyMap_["CN"] = 0.0801;
87 frequencyMap_["NC"] = 0.00521;
88 frequencyMap_["RCHar3"] = -0.00182;
89 frequencyMap_["SigmaN"] = 0.00157;
90 frequencyMap_["PiN"] = -0.00167;
91 frequencyMap_["PiC"] = -0.00896;
92
93 ForceField* forceField_ = info_->getForceField();
94 AtomTypeSet atypes = info_->getSimulatedAtomTypes();
95 PairList* excludes = info_->getExcludedInteractions();
96 int nAtoms =
98
99 RealType rcut;
100 if (info_->getSimParams()->haveCutoffRadius()) {
101 rcut = info_->getSimParams()->getCutoffRadius();
102 } else {
103 rcut = 12.0;
104 }
105
106 EF_ = V3Zero;
107
108 std::vector<RealType> ef;
109 bool efSpec = false;
110
111 if (info_->getSimParams()->haveElectricField()) {
112 efSpec = true;
113 ef = info_->getSimParams()->getElectricField();
114 }
115 if (info_->getSimParams()->haveUniformField()) {
116 efSpec = true;
117 ef = info_->getSimParams()->getUniformField();
118 }
119 if (efSpec) {
120 if (ef.size() != 3) {
121 snprintf(
122 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
123 "NitrileFrequencyMap: Incorrect number of parameters specified for "
124 "uniformField.\n"
125 "\tthere should be 3 parameters, but %lu were specified.\n",
126 ef.size());
127 painCave.isFatal = 1;
128 simError();
129 }
130 EF_.x() = ef[0];
131 EF_.y() = ef[1];
132 EF_.z() = ef[2];
133 }
134
135 excludesForAtom.clear();
136 excludesForAtom.resize(nAtoms);
137
138 for (int i = 0; i < nAtoms; i++) {
139 for (int j = 0; j < nAtoms; j++) {
140 if (excludes->hasPair(i, j)) excludesForAtom[i].push_back(j);
141 }
142 }
143
144 electrostatic_ = new Electrostatic();
145 electrostatic_->setSimInfo(info_);
146 electrostatic_->setForceField(forceField_);
147 electrostatic_->setSimulatedAtomTypes(atypes);
148 electrostatic_->setCutoffRadius(rcut);
149 }
150
151 bool NitrileFrequencyMap::excludeAtomPair(int atom1, int atom2) {
152 for (vector<int>::iterator i = excludesForAtom[atom1].begin();
153 i != excludesForAtom[atom1].end(); ++i) {
154 if ((*i) == atom2) return true;
155 }
156
157 return false;
158 }
159
160 void NitrileFrequencyMap::process() {
161 Molecule* mol;
162 Atom* atom;
163 AtomType* atype;
164 SimInfo::MoleculeIterator mi;
165 Molecule::AtomIterator ai2;
166 Atom* atom2;
167 StuntDouble* sd1;
168 int ii, sdID, molID, sdID2;
169 RealType li(0.0);
170 RealType sPot, s1, s2;
171 RealType freqShift;
172 std::string name;
173 map<string, RealType>::iterator fi;
174 bool excluded;
175 const RealType chrgToKcal = 23.0609;
176
177 DumpReader reader(info_, dumpFilename_);
178 int nFrames = reader.getNFrames();
179
180 nProcessed_ = nFrames / step_;
181
182 std::fill(histogram_.begin(), histogram_.end(), 0.0);
183 std::fill(count_.begin(), count_.end(), 0);
184
185 for (int istep = 0; istep < nFrames; istep += step_) {
186 reader.readFrame(istep);
187 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
188
189 std::fill(freqs_.begin(), freqs_.end(), 0.0);
190
191 if (evaluator1_.isDynamic()) {
192 seleMan1_.setSelectionSet(evaluator1_.evaluate());
193 }
194
195 for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
196 sd1 = seleMan1_.nextSelected(ii)) {
197 sdID = sd1->getGlobalIndex();
198 molID = info_->getGlobalMolMembership(sdID);
199 mol = info_->getMoleculeByGlobalIndex(molID);
200
201 Vector3d CNcentroid = mol->getRigidBodyAt(2)->getPos();
202 Vector3d ra = sd1->getPos();
203
204 atom = dynamic_cast<Atom*>(sd1);
205 atype = atom->getAtomType();
206 name = atype->getName();
207 fi = frequencyMap_.find(name);
208 if (fi != frequencyMap_.end()) {
209 li = (*fi).second;
210 } else {
211 // throw error
212 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
213 "NitrileFrequencyMap::process: Unknown atype requested.\n"
214 "\t(Selection specified %s .)\n",
215 name.c_str());
216 painCave.isFatal = 1;
217 simError();
218 }
219
220 sPot = sd1->getSitePotential();
221
222 // Subtract out the contribution from every other site on this
223 // molecule:
224 for (atom2 = mol->beginAtom(ai2); atom2 != NULL;
225 atom2 = mol->nextAtom(ai2)) {
226 sdID2 = atom2->getGlobalIndex();
227 if (sdID == sdID2) {
228 excluded = true;
229 } else {
230 excluded = excludeAtomPair(sdID, sdID2);
231 }
232
233 electrostatic_->getSitePotentials(atom, atom2, excluded, s1, s2);
234
235 sPot -= s1;
236 }
237
238 // Add the contribution from the electric field:
239
240 sPot += dot(EF_, ra - CNcentroid) * chrgToKcal;
241
242 freqShift = sPot * li;
243
244 // convert the kcal/mol energies to wavenumbers:
245 freqShift *= 349.757;
246
247 freqs_[molID] += freqShift;
248 }
249
250 for (int i = 0; i < info_->getNGlobalMolecules(); ++i) {
251 int binNo =
252 int(nBins_ * (freqs_[i] - minFreq_) / (maxFreq_ - minFreq_));
253
254 count_[binNo]++;
255 }
256 }
257
258 processHistogram();
259 writeProbs();
260 }
261
262 void NitrileFrequencyMap::processHistogram() {
263 int atot = 0;
264 for (unsigned int i = 0; i < count_.size(); ++i)
265 atot += count_[i];
266
267 for (unsigned int i = 0; i < count_.size(); ++i) {
268 histogram_[i] = double(count_[i] / double(atot));
269 }
270 }
271
272 void NitrileFrequencyMap::writeProbs() {
273 std::ofstream rdfStream(outputFilename_.c_str());
274 if (rdfStream.is_open()) {
275 rdfStream << "#NitrileFrequencyMap\n";
276 rdfStream << "#nFrames:\t" << nProcessed_ << "\n";
277 rdfStream << "#selection1: (" << selectionScript1_ << ")";
278 rdfStream << "\n";
279 rdfStream << "#nu\tp(nu))\n";
280 for (unsigned int i = 0; i < histogram_.size(); ++i) {
281 RealType freq = minFreq_ + (RealType)(i) * (maxFreq_ - minFreq_) /
282 (RealType)histogram_.size();
283 rdfStream << freq << "\t" << histogram_[i] << "\n";
284 }
285
286 } else {
287 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
288 "NitrileFrequencyMap: unable to open %s\n",
289 outputFilename_.c_str());
290 painCave.isFatal = 1;
291 simError();
292 }
293
294 rdfStream.close();
295 }
296
297} // namespace OpenMD
AtomType * getAtomType()
Returns the AtomType of this Atom.
Definition A.hpp:93
AtomType is what OpenMD looks to for unchanging data about an atom.
Definition AtomType.hpp:66
PairList class maintains a general purpose list of atom pairs using the global indices of the atoms.
Definition PairList.hpp:63
bool hasPair(int i, int j)
Checks whether pair (i, j) is in this PairList class.
Definition PairList.cpp:140
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
Molecule * getMoleculeByGlobalIndex(int index)
Finds a molecule with a specified global index.
Definition SimInfo.hpp:300
int getNGlobalMolecules()
Returns the total number of molecules in the system.
Definition SimInfo.hpp:126
SnapshotManager * getSnapshotManager()
Returns the snapshot manager.
Definition SimInfo.hpp:248
int getNumberOfAtoms()
Returns the number of atoms.
Definition Snapshot.cpp:202
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.
RealType getSitePotential()
Returns the current site potential of this stuntDouble.
int getGlobalIndex()
Returns the global index 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)