OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
TetrahedralityHBMatrix.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/TetrahedralityHBMatrix.hpp"
46
47#include <vector>
48
49#include "io/DumpReader.hpp"
51#include "utils/Constants.hpp"
52#include "utils/Revision.hpp"
53#include "utils/simError.h"
54
55namespace OpenMD {
56
57 TetrahedralityHBMatrix::TetrahedralityHBMatrix(
58 SimInfo* info, const std::string& filename, const std::string& sele,
59 double rCut, double OOcut, double thetaCut, double OHcut, int nbins) :
60 StaticAnalyser(info, filename, nbins),
61 selectionScript_(sele), seleMan_(info), evaluator_(info), rCut_(rCut),
62 OOCut_(OOcut), thetaCut_(thetaCut), OHCut_(OHcut) {
63 setAnalysisType("Tetrahedrality HBond Matrix");
64 setOutputName(getPrefix(filename) + ".hbq");
65
66 evaluator_.loadScriptString(sele);
67 if (!evaluator_.isDynamic()) {
68 seleMan_.setSelectionSet(evaluator_.evaluate());
69 }
70
71 Q_histogram_.resize(nBins_);
72
73 for (unsigned int i = 0; i < nBins_; i++) {
74 Q_histogram_[i].resize(nBins_);
75 }
76
77 // Q can take values from 0 to 1
78
79 MinQ_ = 0.0;
80 MaxQ_ = 1.1;
81 deltaQ_ = (MaxQ_ - MinQ_) / nBins_;
82
83 std::stringstream params;
84 params << " rCut = " << rCut_ << ", OOcut = " << OOCut_
85 << ", thetacut = " << thetaCut_ << ", OHcut = " << OHCut_
86 << ", nbins = " << nBins_ << ", deltaQ = " << deltaQ_;
87 const std::string paramString = params.str();
88 setParameterString(paramString);
89 }
90
91 void TetrahedralityHBMatrix::initializeHistogram() {
92 for (unsigned int i = 0; i < nBins_; i++) {
93 std::fill(Q_histogram_[i].begin(), Q_histogram_[i].end(), 0);
94 }
95 }
96
97 void TetrahedralityHBMatrix::process() {
98 Molecule* mol1;
99 Molecule* mol2;
100 Molecule* moli;
101 Molecule* molj;
102 SimInfo::MoleculeIterator mi;
103 Vector3d vec;
104 Vector3d ri, rj, rk, rik, rkj, dposition, tposition;
105 RealType r, cospsi;
106 std::vector<Molecule::HBondDonor*>::iterator hbdi;
107 Molecule::HBondDonor* hbd;
108 std::vector<Atom*>::iterator hbai;
109 Atom* hba;
110 Vector3d dPos;
111 Vector3d aPos;
112 Vector3d hPos;
113 Vector3d DH;
114 Vector3d DA;
115 Vector3d HA;
116 Vector3d uDA;
117 RealType DAdist, DHdist, HAdist, theta, ctheta, Qk, q1, q2;
118 std::vector<std::pair<RealType, Molecule*>> myNeighbors;
119 int myIndex, ii, jj, index1, index2;
120
121 bool usePeriodicBoundaryConditions_ =
122 info_->getSimParams()->getUsePeriodicBoundaryConditions();
123
124 DumpReader reader(info_, dumpFilename_);
125 int nFrames = reader.getNFrames();
126
127 for (int istep = 0; istep < nFrames; istep += step_) {
128 reader.readFrame(istep);
129 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
130
131 Q_.clear();
132 Q_.resize(info_->getNGlobalMolecules(), 0.0);
133
134 if (evaluator_.isDynamic()) {
135 seleMan_.setSelectionSet(evaluator_.evaluate());
136 }
137
138 // outer loop is over the selected Molecules:
139
140 for (mol1 = seleMan_.beginSelectedMolecule(ii); mol1 != NULL;
141 mol1 = seleMan_.nextSelectedMolecule(ii)) {
142 myIndex = mol1->getGlobalIndex();
143
144 Qk = 1.0;
145
146 myNeighbors.clear();
147
148 // inner loop is over all Molecules in the system:
149
150 for (mol2 = info_->beginMolecule(mi); mol2 != NULL;
151 mol2 = info_->nextMolecule(mi)) {
152 if (mol2->getGlobalIndex() != myIndex) {
153 vec = mol1->getCom() - mol2->getCom();
154
155 if (usePeriodicBoundaryConditions_)
156 currentSnapshot_->wrapVector(vec);
157
158 r = vec.length();
159
160 // Check to see if neighbor is in bond cutoff
161
162 if (r < rCut_) { myNeighbors.push_back(std::make_pair(r, mol2)); }
163 }
164 }
165
166 // Sort the vector using predicate and std::sort
167 std::sort(myNeighbors.begin(), myNeighbors.end());
168
169 // Use only the 4 closest neighbors to do the rest of the work:
170 // int nbors = myNeighbors.size()> 4 ? 4 : myNeighbors.size();
171 int nbors = myNeighbors.size();
172 int nang = int(0.5 * (nbors * (nbors - 1)));
173
174 rk = mol1->getCom();
175
176 for (int i = 0; i < nbors - 1; i++) {
177 moli = myNeighbors[i].second;
178 ri = moli->getCom();
179 rik = rk - ri;
180 if (usePeriodicBoundaryConditions_) currentSnapshot_->wrapVector(rik);
181
182 rik.normalize();
183
184 for (int j = i + 1; j < nbors; j++) {
185 molj = myNeighbors[j].second;
186 rj = molj->getCom();
187 rkj = rk - rj;
188 if (usePeriodicBoundaryConditions_)
189 currentSnapshot_->wrapVector(rkj);
190 rkj.normalize();
191
192 cospsi = dot(rik, rkj);
193
194 // Calculates scaled Qk for each molecule using calculated
195 // angles from the actual number of nearest neighbors.
196 Qk = Qk - (pow(cospsi + 1.0 / 3.0, 2) * 9.0 / (4.0 * nang));
197 }
198 }
199
200 if (nang > 0) { Q_[myIndex] = Qk; }
201 }
202 // Now to scan for histogram:
203
204 for (mol1 = seleMan_.beginSelectedMolecule(ii); mol1 != NULL;
205 mol1 = seleMan_.nextSelectedMolecule(ii)) {
206 for (mol2 = seleMan_.beginSelectedMolecule(jj); mol2 != NULL;
207 mol2 = seleMan_.nextSelectedMolecule(jj)) {
208 // loop over the possible donors in molecule 1:
209 for (hbd = mol1->beginHBondDonor(hbdi); hbd != NULL;
210 hbd = mol1->nextHBondDonor(hbdi)) {
211 dPos = hbd->donorAtom->getPos();
212 hPos = hbd->donatedHydrogen->getPos();
213 DH = hPos - dPos;
214 currentSnapshot_->wrapVector(DH);
215 DHdist = DH.length();
216
217 // loop over the possible acceptors in molecule 2:
218 for (hba = mol2->beginHBondAcceptor(hbai); hba != NULL;
219 hba = mol2->nextHBondAcceptor(hbai)) {
220 aPos = hba->getPos();
221 DA = aPos - dPos;
222 currentSnapshot_->wrapVector(DA);
223 DAdist = DA.length();
224
225 // Distance criteria: are the donor and acceptor atoms
226 // close enough?
227 if (DAdist < OOCut_) {
228 HA = aPos - hPos;
229 currentSnapshot_->wrapVector(HA);
230 HAdist = HA.length();
231
232 ctheta = dot(DH, DA) / (DHdist * DAdist);
233 theta = acos(ctheta) * 180.0 / Constants::PI;
234
235 // Angle criteria: are the D-H and D-A and vectors close?
236 if (theta < thetaCut_ && HAdist < OHCut_) {
237 // We have a hydrogen bond!
238 index1 = mol1->getGlobalIndex();
239 index2 = mol2->getGlobalIndex();
240 q1 = Q_[index1];
241 q2 = Q_[index2];
242 collectHistogram(q1, q2);
243 }
244 }
245 }
246 }
247 // now loop over the possible acceptors in molecule 1:
248 for (hba = mol1->beginHBondAcceptor(hbai); hba != NULL;
249 hba = mol1->nextHBondAcceptor(hbai)) {
250 aPos = hba->getPos();
251
252 // loop over the possible donors in molecule 2:
253 for (hbd = mol2->beginHBondDonor(hbdi); hbd != NULL;
254 hbd = mol2->nextHBondDonor(hbdi)) {
255 dPos = hbd->donorAtom->getPos();
256
257 DA = aPos - dPos;
258 currentSnapshot_->wrapVector(DA);
259 DAdist = DA.length();
260
261 // Distance criteria: are the donor and acceptor atoms
262 // close enough?
263 if (DAdist < OOCut_) {
264 hPos = hbd->donatedHydrogen->getPos();
265 HA = aPos - hPos;
266 currentSnapshot_->wrapVector(HA);
267 HAdist = HA.length();
268
269 DH = hPos - dPos;
270 currentSnapshot_->wrapVector(DH);
271 DHdist = DH.length();
272 ctheta = dot(DH, DA) / (DHdist * DAdist);
273 theta = acos(ctheta) * 180.0 / Constants::PI;
274 // Angle criteria: are the D-H and D-A and vectors close?
275 if (theta < thetaCut_ && HAdist < OHCut_) {
276 // We have a hydrogen bond!
277 // keep these indexed by donor then acceptor:
278 index1 = mol2->getGlobalIndex();
279 index2 = mol1->getGlobalIndex();
280 q1 = Q_[index1];
281 q2 = Q_[index2];
282 collectHistogram(q1, q2);
283 }
284 }
285 }
286 }
287 }
288 }
289 }
290 writeOutput();
291 }
292
293 void TetrahedralityHBMatrix::collectHistogram(RealType q1, RealType q2) {
294 if (q1 > MinQ_ && q1 < MaxQ_) {
295 int bin1 = int((q1 - MinQ_) / deltaQ_);
296 if (q2 > MinQ_ && q2 < MaxQ_) {
297 int bin2 = int((q2 - MinQ_) / deltaQ_);
298 Q_histogram_[bin1][bin2] += 1;
299 count_++;
300 } else {
301 // cerr << "q2 = " << q2 << "\n";
302 }
303 } else {
304 // cerr << "q1 = " << q1 << "\n";
305 }
306 }
307
308 void TetrahedralityHBMatrix::writeOutput() {
309 std::ofstream ofs(outputFilename_.c_str());
310 if (ofs.is_open()) {
311 Revision r;
312 ofs << "# " << getAnalysisType() << "\n";
313 ofs << "# OpenMD " << r.getFullRevision() << "\n";
314 ofs << "# " << r.getBuildDate() << "\n";
315 ofs << "# selection script: \"" << selectionScript_ << "\"\n";
316
317 if (!paramString_.empty())
318 ofs << "# parameters: " << paramString_ << "\n";
319
320 ofs << "# Hydrogen Bond count: " << count_ << "\n";
321
322 for (unsigned int i = 0; i < Q_histogram_.size(); ++i) {
323 for (unsigned int j = 0; j < Q_histogram_[i].size(); ++j) {
324 ofs << RealType(Q_histogram_[i][j]) / RealType(count_) << "\t";
325 }
326 ofs << "\n";
327 }
328
329 } else {
330 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
331 "TetrahedralityHBMatrix: unable to open %s\n",
332 outputFilename_.c_str());
333 painCave.isFatal = 1;
334 simError();
335 }
336
337 ofs.close();
338 }
339} // 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)