OpenMD 3.1
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
HBondJump.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/dynamicProps/HBondJump.hpp"
46
47#include <algorithm>
48
49#include "utils/Constants.hpp"
50#include "utils/Revision.hpp"
51
52namespace OpenMD {
53
54 HBondJump::HBondJump(SimInfo* info, const std::string& filename,
55 const std::string& sele1, const std::string& sele2,
56 double OOcut, double thetaCut, double OHcut) :
57 TimeCorrFunc<RealType>(info, filename, sele1, sele2),
58 OOCut_(OOcut), thetaCut_(thetaCut), OHCut_(OHcut),
59 sele1_minus_common_(info), sele2_minus_common_(info), common_(info) {
60 setCorrFuncType("HBondJump");
61 setOutputName(getPrefix(dumpFilename_) + ".jump");
62
63 std::stringstream params;
64 params << " OOcut = " << OOCut_ << ", thetacut = " << thetaCut_
65 << ", OHcut = " << OHCut_;
66 const std::string paramString = params.str();
67 setParameterString(paramString);
68
69 if (!uniqueSelections_) { seleMan2_ = seleMan1_; }
70 evaluator1_.loadScriptString(selectionScript1_);
71 // if selection is static, we only need to evaluate it once
72 if (!evaluator1_.isDynamic()) {
73 seleMan1_.setSelectionSet(evaluator1_.evaluate());
74 validateSelection(seleMan1_);
75 }
76
77 if (uniqueSelections_) {
78 evaluator2_.loadScriptString(selectionScript2_);
79 if (!evaluator2_.isDynamic()) {
80 seleMan2_.setSelectionSet(evaluator2_.evaluate());
81 validateSelection(seleMan2_);
82 }
83 }
84
85 if (!evaluator1_.isDynamic() && !evaluator2_.isDynamic()) {
86 // If all selections are static, we can pre-set the selections:
87 common_ = seleMan1_ & seleMan2_;
88 sele1_minus_common_ = seleMan1_ - common_;
89 sele2_minus_common_ = seleMan2_ - common_;
90 }
91
92 // nFrames_ is initialized in MultipassCorrFunc:
93 GIDtoH_.resize(nFrames_);
94 hydrogen_.resize(nFrames_);
95 acceptor_.resize(nFrames_);
96 lastAcceptor_.resize(nFrames_);
97 acceptorStartFrame_.resize(nFrames_);
98 selected_.resize(nFrames_);
99 }
100
101 void HBondJump::computeFrame(int istep) {
102 // Map of atomic global IDs to HBond donor hydrogens:
103 GIDtoH_[istep].resize(info_->getNGlobalAtoms(), -1);
104
105 // Find all of the HBonds in this frame
106 findHBonds(istep);
107
108 if (!uniqueSelections_) { seleMan2_ = seleMan1_; }
109 if (evaluator1_.isDynamic()) {
110 seleMan1_.setSelectionSet(evaluator1_.evaluate());
111 }
112 if (uniqueSelections_ && evaluator2_.isDynamic()) {
113 seleMan2_.setSelectionSet(evaluator2_.evaluate());
114 }
115 if (evaluator1_.isDynamic() || evaluator2_.isDynamic()) {
116 common_ = seleMan1_ & seleMan2_;
117 sele1_minus_common_ = seleMan1_ - common_;
118 sele2_minus_common_ = seleMan2_ - common_;
119 }
120
121 // Label the found HBonds as selected:
122 processNonOverlapping(istep, sele1_minus_common_, seleMan2_);
123 processNonOverlapping(istep, common_, sele2_minus_common_);
124 processOverlapping(istep, common_);
125 }
126
127 void HBondJump::correlation() {
128 std::vector<int> s1;
129 std::vector<int>::iterator i1;
130 RealType corrVal;
131 int index1, index2, count, gid, aInd1, aInd2;
132
133 for (int i = 0; i < nFrames_; ++i) {
134 RealType time1 = times_[i];
135 s1 = hydrogen_[i];
136
137 for (int j = i; j < nFrames_; ++j) {
138 // Perform a sanity check on the actual configuration times to
139 // make sure the configurations are spaced the same amount the
140 // sample time said they were spaced:
141
142 RealType time2 = times_[j];
143
144 if (fabs((time2 - time1) - (j - i) * deltaTime_) > 1.0e-4) {
145 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
146 "HBondJump::correlation Error: sampleTime (%f)\n"
147 "\tin %s does not match actual time-spacing between\n"
148 "\tconfigurations %d (t = %f) and %d (t = %f).\n",
149 deltaTime_, dumpFilename_.c_str(), i, time1, j, time2);
150 painCave.isFatal = 1;
151 simError();
152 }
153
154 int timeBin = int((time2 - time1) / deltaTime_ + 0.5);
155
156 corrVal = 0.0;
157 count = 0;
158
159 // loop over the Hydrogens found in frame i:
160
161 for (i1 = s1.begin(); i1 != s1.end(); ++i1) {
162 // gid is the global ID of Hydrogen index1 in frame i:
163 gid = *i1;
164 index1 = GIDtoH_[i][gid];
165
166 // find matching hydrogen in frame j:
167 index2 = GIDtoH_[j][gid];
168
169 if (selected_[i][index1]) {
170 count++;
171
172 if (acceptor_[i][index1] == -1) {
173 aInd1 = lastAcceptor_[i][index1];
174 } else {
175 aInd1 = acceptor_[i][index1];
176 }
177
178 if (acceptor_[j][index2] == -1) {
179 aInd2 = lastAcceptor_[j][index2];
180 } else {
181 aInd2 = acceptor_[j][index2];
182 }
183
184 // aInd1 = acceptor_[i][index1];
185 // aInd2 = acceptor_[j][index2];
186
187 if (aInd1 != aInd2) {
188 // different acceptor so nA(0) . nB(t) = 1
189 corrVal += 1;
190 } else {
191 // same acceptor, but we need to look at the start frames
192 // for these H-bonds to make sure it is the same H-bond:
193 if (acceptorStartFrame_[i][index1] !=
194 acceptorStartFrame_[j][index2]) {
195 // different start frame, so this is considered a
196 // different H-bond:
197 corrVal += 1;
198 } else {
199 // same start frame, so this is considered the same H-bond:
200 corrVal += 0;
201 }
202 }
203 }
204 }
205 histogram_[timeBin] += corrVal;
206 count_[timeBin] += count;
207 }
208 }
209 }
210
211 void HBondJump::postCorrelate() {
212 for (unsigned int i = 0; i < nTimeBins_; ++i) {
213 if (count_[i] > 0) {
214 histogram_[i] /= count_[i];
215 } else {
216 histogram_[i] = 0;
217 }
218 histogram_[i] = 1.0 - histogram_[i];
219 }
220 }
221
222 void HBondJump::processNonOverlapping(int frame, SelectionManager& sman1,
223 SelectionManager& sman2) {
224 Molecule* mol1;
225 Molecule* mol2;
226 int i, j;
227 std::vector<Molecule::HBondDonor*>::iterator hbdi;
228 Molecule::HBondDonor* hbd;
229 std::vector<Atom*>::iterator hbai;
230 Atom* hba;
231 int hInd, index, aInd;
232
233 // This is the same as a non-overlapping pairwise loop structure:
234 // for (int i = 0; i < ni ; ++i ) {
235 // for (int j = 0; j < nj; ++j) {}
236 // }
237
238 for (mol1 = sman1.beginSelectedMolecule(i); mol1 != NULL;
239 mol1 = sman1.nextSelectedMolecule(i)) {
240 for (mol2 = sman2.beginSelectedMolecule(j); mol2 != NULL;
241 mol2 = sman2.nextSelectedMolecule(j)) {
242 // loop over the possible donors in molecule 1:
243 for (hbd = mol1->beginHBondDonor(hbdi); hbd != NULL;
244 hbd = mol1->nextHBondDonor(hbdi)) {
245 hInd = hbd->donatedHydrogen->getGlobalIndex();
246 index = GIDtoH_[frame][hInd];
247 aInd = acceptor_[frame][index];
248
249 for (hba = mol2->beginHBondAcceptor(hbai); hba != NULL;
250 hba = mol2->nextHBondAcceptor(hbai)) {
251 if (hba->getGlobalIndex() == aInd) {
252 selected_[frame][index] = true;
253 }
254 }
255 }
256
257 // loop over the possible donors in molecule 2:
258 for (hbd = mol2->beginHBondDonor(hbdi); hbd != NULL;
259 hbd = mol2->nextHBondDonor(hbdi)) {
260 hInd = hbd->donatedHydrogen->getGlobalIndex();
261 index = GIDtoH_[frame][hInd];
262 aInd = acceptor_[frame][index];
263
264 for (hba = mol1->beginHBondAcceptor(hbai); hba != NULL;
265 hba = mol1->nextHBondAcceptor(hbai)) {
266 if (hba->getGlobalIndex() == aInd) {
267 selected_[frame][index] = true;
268 }
269 }
270 }
271 }
272 }
273 }
274
275 void HBondJump::processOverlapping(int frame, SelectionManager& sman) {
276 Molecule* mol1;
277 Molecule* mol2;
278 int i, j;
279
280 std::vector<Molecule::HBondDonor*>::iterator hbdi;
281 Molecule::HBondDonor* hbd;
282 std::vector<Atom*>::iterator hbai;
283 Atom* hba;
284 int hInd, index, aInd;
285
286 // This is the same as a pairwise loop structure:
287 // for (int i = 0; i < n-1 ; ++i ) {
288 // for (int j = i + 1; j < n; ++j) {}
289 // }
290
291 for (mol1 = sman.beginSelectedMolecule(i); mol1 != NULL;
292 mol1 = sman.nextSelectedMolecule(i)) {
293 // loop over the possible donors in molecule 1:
294 for (hbd = mol1->beginHBondDonor(hbdi); hbd != NULL;
295 hbd = mol1->nextHBondDonor(hbdi)) {
296 hInd = hbd->donatedHydrogen->getGlobalIndex();
297 index = GIDtoH_[frame][hInd];
298 aInd = acceptor_[frame][index];
299
300 for (j = i, mol2 = sman.nextSelectedMolecule(j); mol2 != NULL;
301 mol2 = sman.nextSelectedMolecule(j)) {
302 for (hba = mol2->beginHBondAcceptor(hbai); hba != NULL;
303 hba = mol2->nextHBondAcceptor(hbai)) {
304 if (hba->getGlobalIndex() == aInd) {
305 selected_[frame][index] = true;
306 }
307 }
308 }
309 }
310
311 for (hba = mol1->beginHBondAcceptor(hbai); hba != NULL;
312 hba = mol1->nextHBondAcceptor(hbai)) {
313 aInd = hba->getGlobalIndex();
314
315 for (j = i, mol2 = sman.nextSelectedMolecule(j); mol2 != NULL;
316 mol2 = sman.nextSelectedMolecule(j)) {
317 for (hbd = mol2->beginHBondDonor(hbdi); hbd != NULL;
318 hbd = mol2->nextHBondDonor(hbdi)) {
319 hInd = hbd->donatedHydrogen->getGlobalIndex();
320 index = GIDtoH_[frame][hInd];
321
322 if (acceptor_[frame][index] == aInd) {
323 selected_[frame][index] = true;
324 }
325 }
326 }
327 }
328 }
329 }
330
331 void HBondJump::findHBonds(int frame) {
332 Molecule* mol1;
333 Molecule* mol2;
334 SimInfo::MoleculeIterator mi, mj;
335 std::vector<Molecule::HBondDonor*>::iterator hbdi;
336 Molecule::HBondDonor* hbd;
337 std::vector<Atom*>::iterator hbai;
338 Atom* hba;
339 Vector3d dPos, hPos, aPos;
340 int hInd, index, aInd;
341
342 // Register all the possible HBond donor hydrogens:
343 for (mol1 = info_->beginMolecule(mi); mol1 != NULL;
344 mol1 = info_->nextMolecule(mi)) {
345 for (hbd = mol1->beginHBondDonor(hbdi); hbd != NULL;
346 hbd = mol1->nextHBondDonor(hbdi)) {
347 hInd = hbd->donatedHydrogen->getGlobalIndex();
348 index = registerHydrogen(frame, hInd);
349 }
350 }
351
352 for (mol1 = info_->beginMolecule(mi); mol1 != NULL;
353 mol1 = info_->nextMolecule(mi)) {
354 for (hbd = mol1->beginHBondDonor(hbdi); hbd != NULL;
355 hbd = mol1->nextHBondDonor(hbdi)) {
356 hInd = hbd->donatedHydrogen->getGlobalIndex();
357 index = GIDtoH_[frame][hInd];
358
359 dPos = hbd->donorAtom->getPos();
360 hPos = hbd->donatedHydrogen->getPos();
361
362 for (mj = mi, mol2 = info_->beginMolecule(mj); mol2 != NULL;
363 mol2 = info_->nextMolecule(mj)) {
364 for (hba = mol2->beginHBondAcceptor(hbai); hba != NULL;
365 hba = mol2->nextHBondAcceptor(hbai)) {
366 aPos = hba->getPos();
367
368 if (isHBond(dPos, hPos, aPos)) {
369 aInd = hba->getGlobalIndex();
370 registerHydrogenBond(frame, index, hInd, aInd);
371 }
372 }
373 }
374 }
375
376 for (hba = mol1->beginHBondAcceptor(hbai); hba != NULL;
377 hba = mol1->nextHBondAcceptor(hbai)) {
378 aPos = hba->getPos();
379
380 for (mj = mi, mol2 = info_->beginMolecule(mj); mol2 != NULL;
381 mol2 = info_->nextMolecule(mj)) {
382 for (hbd = mol2->beginHBondDonor(hbdi); hbd != NULL;
383 hbd = mol2->nextHBondDonor(hbdi)) {
384 hInd = hbd->donatedHydrogen->getGlobalIndex();
385 // no need to register, just look up the index:
386 index = GIDtoH_[frame][hInd];
387
388 dPos = hbd->donorAtom->getPos();
389 hPos = hbd->donatedHydrogen->getPos();
390
391 if (isHBond(dPos, hPos, aPos)) {
392 aInd = hba->getGlobalIndex();
393 registerHydrogenBond(frame, index, hInd, aInd);
394 }
395 }
396 }
397 }
398 }
399 }
400
401 bool HBondJump::isHBond(Vector3d donorPos, Vector3d hydrogenPos,
402 Vector3d acceptorPos) {
403 Vector3d DA = acceptorPos - donorPos;
404 currentSnapshot_->wrapVector(DA);
405 RealType DAdist = DA.length();
406
407 // Distance criteria: are the donor and acceptor atoms
408 // close enough?
409 if (DAdist < OOCut_) {
410 Vector3d DH = hydrogenPos - donorPos;
411 currentSnapshot_->wrapVector(DH);
412 RealType DHdist = DH.length();
413
414 Vector3d HA = acceptorPos - hydrogenPos;
415 currentSnapshot_->wrapVector(HA);
416 RealType HAdist = HA.length();
417
418 RealType ctheta = dot(DH, DA) / (DHdist * DAdist);
419 RealType theta = acos(ctheta) * 180.0 / Constants::PI;
420
421 // Angle criteria: are the D-H and D-A and vectors close?
422 if (theta < thetaCut_ && HAdist < OHCut_) { return true; }
423 }
424 return false;
425 }
426
427 int HBondJump::registerHydrogen(int frame, int hIndex) {
428 int index;
429
430 // If this hydrogen wasn't already registered, register it:
431 if (GIDtoH_[frame][hIndex] == -1) {
432 index = hydrogen_[frame].size();
433 GIDtoH_[frame][hIndex] = index;
434 hydrogen_[frame].push_back(hIndex);
435 acceptor_[frame].push_back(-1);
436 selected_[frame].push_back(false);
437
438 if (frame == 0) {
439 lastAcceptor_[frame].push_back(-1);
440 acceptorStartFrame_[frame].push_back(frame);
441 } else {
442 // Copy the last acceptor.
443 int prevIndex = GIDtoH_[frame - 1][hIndex];
444 lastAcceptor_[frame].push_back(lastAcceptor_[frame - 1][prevIndex]);
445 acceptorStartFrame_[frame].push_back(
446 acceptorStartFrame_[frame - 1][prevIndex]);
447 }
448 } else {
449 // This hydrogen was already registered. Just return the index:
450 index = GIDtoH_[frame][hIndex];
451 }
452 return index;
453 }
454
455 void HBondJump::registerHydrogenBond(int frame, int index, int hIndex,
456 int acceptorIndex) {
457 acceptor_[frame][index] = acceptorIndex;
458 lastAcceptor_[frame][index] = acceptorIndex;
459
460 if (frame == 0) {
461 acceptorStartFrame_[frame][index] = frame;
462 } else {
463 int prevIndex = GIDtoH_[frame - 1][hIndex];
464 if (acceptorIndex != lastAcceptor_[frame - 1][prevIndex]) {
465 acceptorStartFrame_[frame][index] = frame;
466 } else {
467 acceptorStartFrame_[frame][index] =
468 acceptorStartFrame_[frame - 1][prevIndex];
469 }
470 }
471 }
472
473 HBondJumpZ::HBondJumpZ(SimInfo* info, const std::string& filename,
474 const std::string& sele1, const std::string& sele2,
475 double OOcut, double thetaCut, double OHcut,
476 int nZBins, int axis) :
477 HBondJump(info, filename, sele1, sele2, OOcut, thetaCut, OHcut),
478 nZBins_(nZBins), axis_(axis) {
479 setCorrFuncType("HBondJumpZ");
480 setOutputName(getPrefix(dumpFilename_) + ".jumpZ");
481
482 switch (axis_) {
483 case 0:
484 axisLabel_ = "x";
485 break;
486 case 1:
487 axisLabel_ = "y";
488 break;
489 case 2:
490 default:
491 axisLabel_ = "z";
492 break;
493 }
494
495 zbin_.resize(nFrames_);
496 histogram_.resize(nTimeBins_);
497 counts_.resize(nTimeBins_);
498 for (unsigned int i = 0; i < nTimeBins_; i++) {
499 histogram_[i].resize(nZBins_);
500 std::fill(histogram_[i].begin(), histogram_[i].end(), 0.0);
501 counts_[i].resize(nZBins_);
502 std::fill(counts_[i].begin(), counts_[i].end(), 0);
503 }
504 }
505
506 void HBondJumpZ::findHBonds(int frame) {
507 Molecule* mol1;
508 Molecule* mol2;
509 SimInfo::MoleculeIterator mi, mj;
510 std::vector<Molecule::HBondDonor*>::iterator hbdi;
511 Molecule::HBondDonor* hbd;
512 std::vector<Atom*>::iterator hbai;
513 Atom* hba;
514 Vector3d dPos, hPos, aPos, pos;
515 int hInd, index, aInd, zBin;
516
517 Mat3x3d hmat = currentSnapshot_->getHmat();
518 RealType halfBoxZ_ = hmat(axis_, axis_) / 2.0;
519
520 // Register all the possible HBond donor hydrogens:
521 for (mol1 = info_->beginMolecule(mi); mol1 != NULL;
522 mol1 = info_->nextMolecule(mi)) {
523 for (hbd = mol1->beginHBondDonor(hbdi); hbd != NULL;
524 hbd = mol1->nextHBondDonor(hbdi)) {
525 hInd = hbd->donatedHydrogen->getGlobalIndex();
526 index = registerHydrogen(frame, hInd);
527 }
528 }
529
530 for (mol1 = info_->beginMolecule(mi); mol1 != NULL;
531 mol1 = info_->nextMolecule(mi)) {
532 for (hbd = mol1->beginHBondDonor(hbdi); hbd != NULL;
533 hbd = mol1->nextHBondDonor(hbdi)) {
534 hInd = hbd->donatedHydrogen->getGlobalIndex();
535 index = GIDtoH_[frame][hInd];
536
537 dPos = hbd->donorAtom->getPos();
538 hPos = hbd->donatedHydrogen->getPos();
539
540 for (mj = mi, mol2 = info_->beginMolecule(mj); mol2 != NULL;
541 mol2 = info_->nextMolecule(mj)) {
542 for (hba = mol2->beginHBondAcceptor(hbai); hba != NULL;
543 hba = mol2->nextHBondAcceptor(hbai)) {
544 aPos = hba->getPos();
545
546 if (isHBond(dPos, hPos, aPos)) {
547 aInd = hba->getGlobalIndex();
548 registerHydrogenBond(frame, index, hInd, aInd);
549 pos = hPos;
550 if (info_->getSimParams()->getUsePeriodicBoundaryConditions())
551 currentSnapshot_->wrapVector(pos);
552 zBin =
553 int(nZBins_ * (halfBoxZ_ + pos[axis_]) / hmat(axis_, axis_));
554 zbin_[frame][index] = zBin;
555 }
556 }
557 }
558 }
559
560 for (hba = mol1->beginHBondAcceptor(hbai); hba != NULL;
561 hba = mol1->nextHBondAcceptor(hbai)) {
562 aPos = hba->getPos();
563
564 for (mj = mi, mol2 = info_->beginMolecule(mj); mol2 != NULL;
565 mol2 = info_->nextMolecule(mj)) {
566 for (hbd = mol2->beginHBondDonor(hbdi); hbd != NULL;
567 hbd = mol2->nextHBondDonor(hbdi)) {
568 hInd = hbd->donatedHydrogen->getGlobalIndex();
569 // no need to register, just look up the index:
570 index = GIDtoH_[frame][hInd];
571
572 dPos = hbd->donorAtom->getPos();
573 hPos = hbd->donatedHydrogen->getPos();
574
575 if (isHBond(dPos, hPos, aPos)) {
576 aInd = hba->getGlobalIndex();
577 registerHydrogenBond(frame, index, hInd, aInd);
578 pos = hPos;
579 if (info_->getSimParams()->getUsePeriodicBoundaryConditions())
580 currentSnapshot_->wrapVector(pos);
581 zBin =
582 int(nZBins_ * (halfBoxZ_ + pos[axis_]) / hmat(axis_, axis_));
583 zbin_[frame][index] = zBin;
584 }
585 }
586 }
587 }
588 }
589 }
590
591 int HBondJumpZ::registerHydrogen(int frame, int hIndex) {
592 int index;
593
594 // If this hydrogen wasn't already registered, register it:
595 if (GIDtoH_[frame][hIndex] == -1) {
596 index = hydrogen_[frame].size();
597 GIDtoH_[frame][hIndex] = index;
598 hydrogen_[frame].push_back(hIndex);
599 acceptor_[frame].push_back(-1);
600 selected_[frame].push_back(false);
601 zbin_[frame].push_back(-1);
602
603 if (frame == 0) {
604 lastAcceptor_[frame].push_back(-1);
605 acceptorStartFrame_[frame].push_back(frame);
606 } else {
607 // Copy the last acceptor.
608 int prevIndex = GIDtoH_[frame - 1][hIndex];
609 lastAcceptor_[frame].push_back(lastAcceptor_[frame - 1][prevIndex]);
610 acceptorStartFrame_[frame].push_back(
611 acceptorStartFrame_[frame - 1][prevIndex]);
612 }
613 } else {
614 // This hydrogen was already registered. Just return the index:
615 index = GIDtoH_[frame][hIndex];
616 }
617 return index;
618 }
619
620 void HBondJumpZ::correlation() {
621 std::vector<int> s1;
622 std::vector<int>::iterator i1;
623 int index1, index2, gid, aInd1, aInd2, zBin;
624
625 for (int i = 0; i < nFrames_; ++i) {
626 RealType time1 = times_[i];
627 s1 = hydrogen_[i];
628
629 for (int j = i; j < nFrames_; ++j) {
630 // Perform a sanity check on the actual configuration times to
631 // make sure the configurations are spaced the same amount the
632 // sample time said they were spaced:
633
634 RealType time2 = times_[j];
635
636 if (fabs((time2 - time1) - (j - i) * deltaTime_) > 1.0e-4) {
637 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
638 "HBondJump::correlation Error: sampleTime (%f)\n"
639 "\tin %s does not match actual time-spacing between\n"
640 "\tconfigurations %d (t = %f) and %d (t = %f).\n",
641 deltaTime_, dumpFilename_.c_str(), i, time1, j, time2);
642 painCave.isFatal = 1;
643 simError();
644 }
645
646 int timeBin = int((time2 - time1) / deltaTime_ + 0.5);
647
648 // loop over the Hydrogens found in frame i:
649
650 for (i1 = s1.begin(); i1 != s1.end(); ++i1) {
651 // gid is the global ID of Hydrogen index1 in frame i:
652 gid = *i1;
653 index1 = GIDtoH_[i][gid];
654
655 // find matching hydrogen in frame j:
656 index2 = GIDtoH_[j][gid];
657
658 if (selected_[i][index1]) {
659 zBin = zbin_[i][index1];
660 counts_[timeBin][zBin]++;
661
662 if (acceptor_[i][index1] == -1) {
663 aInd1 = lastAcceptor_[i][index1];
664 } else {
665 aInd1 = acceptor_[i][index1];
666 }
667
668 if (acceptor_[j][index2] == -1) {
669 aInd2 = lastAcceptor_[j][index2];
670 } else {
671 aInd2 = acceptor_[j][index2];
672 }
673
674 // aInd1 = acceptor_[i][index1];
675 // aInd2 = acceptor_[j][index2];
676
677 if (aInd1 != aInd2) {
678 // different acceptor so nA(0) . nB(t) = 1
679 histogram_[timeBin][zBin] += 1;
680 } else {
681 // same acceptor, but we need to look at the start frames
682 // for these H-bonds to make sure it is the same H-bond:
683 if (acceptorStartFrame_[i][index1] !=
684 acceptorStartFrame_[j][index2]) {
685 // different start frame, so this is considered a
686 // different H-bond:
687 histogram_[timeBin][zBin] += 1;
688 } else {
689 // same start frame, so this is considered the same H-bond:
690 histogram_[timeBin][zBin] += 0;
691 }
692 }
693 }
694 }
695 }
696 }
697 }
698
699 void HBondJumpZ::postCorrelate() {
700 for (unsigned int i = 0; i < nTimeBins_; ++i) {
701 for (unsigned int j = 0; j < nZBins_; ++j) {
702 if (counts_[i][j] > 0) {
703 histogram_[i][j] /= counts_[i][j];
704 } else {
705 histogram_[i][j] = 0;
706 }
707 histogram_[i][j] = 1.0 - histogram_[i][j];
708 }
709 }
710 }
711
712 void HBondJumpZ::writeCorrelate() {
713 ofstream ofs(outputFilename_.c_str());
714
715 if (ofs.is_open()) {
716 Revision r;
717
718 ofs << "# " << getCorrFuncType() << "\n";
719 ofs << "# OpenMD " << r.getFullRevision() << "\n";
720 ofs << "# " << r.getBuildDate() << "\n";
721 ofs << "# selection script1: \"" << selectionScript1_;
722 ofs << "\"\tselection script2: \"" << selectionScript2_ << "\"\n";
723 ofs << "# privilegedAxis computed as " << axisLabel_ << " axis \n";
724 if (!paramString_.empty())
725 ofs << "# parameters: " << paramString_ << "\n";
726 ofs << "#time\tcorrVal\n";
727
728 for (unsigned int i = 0; i < nTimeBins_; ++i) {
729 ofs << times_[i] - times_[0];
730
731 for (unsigned int j = 0; j < nZBins_; ++j) {
732 ofs << "\t" << histogram_[i][j];
733 }
734 ofs << "\n";
735 }
736
737 } else {
738 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
739 "HBondJumpZ::writeCorrelate Error: fail to open %s\n",
740 outputFilename_.c_str());
741 painCave.isFatal = 1;
742 simError();
743 }
744
745 ofs.close();
746 }
747
748 HBondJumpR::HBondJumpR(SimInfo* info, const std::string& filename,
749 const std::string& sele1, const std::string& sele2,
750 const std::string& sele3, RealType OOcut,
751 RealType thetaCut, RealType OHcut, RealType len,
752 int nRBins) :
753 HBondJump(info, filename, sele1, sele2, OOcut, thetaCut, OHcut),
754 len_(len), nRBins_(nRBins), seleMan3_(info_), selectionScript3_(sele3),
755 evaluator3_(info_) {
756 setCorrFuncType("HBondJumpR");
757 setOutputName(getPrefix(dumpFilename_) + ".jumpR");
758
759 rbin_.resize(nFrames_);
760 histogram_.resize(nTimeBins_);
761 counts_.resize(nTimeBins_);
762 for (unsigned int i = 0; i < nTimeBins_; i++) {
763 histogram_[i].resize(nRBins_);
764 counts_[i].resize(nRBins_);
765 }
766
767 evaluator3_.loadScriptString(selectionScript3_);
768 if (!evaluator3_.isDynamic()) {
769 seleMan3_.setSelectionSet(evaluator3_.evaluate());
770 }
771
772 deltaR_ = len_ / nRBins_;
773 }
774
775 void HBondJumpR::findHBonds(int frame) {
776 Molecule* mol1;
777 Molecule* mol2;
778 SimInfo::MoleculeIterator mi, mj;
779 std::vector<Molecule::HBondDonor*>::iterator hbdi;
780 Molecule::HBondDonor* hbd;
781 std::vector<Atom*>::iterator hbai;
782 Atom* hba;
783 Vector3d dPos, hPos, aPos, pos;
784 int hInd, index, aInd;
785 StuntDouble* sd3;
786 Vector3d vec;
787 RealType r;
788 int isd3;
789
790 if (evaluator3_.isDynamic()) {
791 seleMan3_.setSelectionSet(evaluator3_.evaluate());
792 }
793
794 bool usePeriodicBoundaryConditions_ =
795 info_->getSimParams()->getUsePeriodicBoundaryConditions();
796
797 // Register all the possible HBond donor hydrogens:
798 for (mol1 = info_->beginMolecule(mi); mol1 != NULL;
799 mol1 = info_->nextMolecule(mi)) {
800 for (hbd = mol1->beginHBondDonor(hbdi); hbd != NULL;
801 hbd = mol1->nextHBondDonor(hbdi)) {
802 hInd = hbd->donatedHydrogen->getGlobalIndex();
803 index = registerHydrogen(frame, hInd);
804 }
805 }
806
807 for (mol1 = info_->beginMolecule(mi); mol1 != NULL;
808 mol1 = info_->nextMolecule(mi)) {
809 for (hbd = mol1->beginHBondDonor(hbdi); hbd != NULL;
810 hbd = mol1->nextHBondDonor(hbdi)) {
811 hInd = hbd->donatedHydrogen->getGlobalIndex();
812 index = GIDtoH_[frame][hInd];
813
814 dPos = hbd->donorAtom->getPos();
815 hPos = hbd->donatedHydrogen->getPos();
816
817 for (mj = mi, mol2 = info_->beginMolecule(mj); mol2 != NULL;
818 mol2 = info_->nextMolecule(mj)) {
819 for (hba = mol2->beginHBondAcceptor(hbai); hba != NULL;
820 hba = mol2->nextHBondAcceptor(hbai)) {
821 aPos = hba->getPos();
822
823 if (isHBond(dPos, hPos, aPos)) {
824 aInd = hba->getGlobalIndex();
825 registerHydrogenBond(frame, index, hInd, aInd);
826 pos = hPos;
827
828 RealType shortest = HONKING_LARGE_VALUE;
829
830 // loop over selection 3 to find closest atom in selection 3:
831 for (sd3 = seleMan3_.beginSelected(isd3); sd3 != NULL;
832 sd3 = seleMan3_.nextSelected(isd3)) {
833 vec = pos - sd3->getPos();
834
835 if (usePeriodicBoundaryConditions_)
836 currentSnapshot_->wrapVector(vec);
837
838 r = vec.length();
839
840 if (r < shortest) shortest = r;
841 }
842
843 int whichBin = int(shortest / deltaR_);
844 if (whichBin < int(nRBins_)) { rbin_[frame][index] = whichBin; }
845 }
846 }
847 }
848 }
849
850 for (hba = mol1->beginHBondAcceptor(hbai); hba != NULL;
851 hba = mol1->nextHBondAcceptor(hbai)) {
852 aPos = hba->getPos();
853
854 for (mj = mi, mol2 = info_->beginMolecule(mj); mol2 != NULL;
855 mol2 = info_->nextMolecule(mj)) {
856 for (hbd = mol2->beginHBondDonor(hbdi); hbd != NULL;
857 hbd = mol2->nextHBondDonor(hbdi)) {
858 hInd = hbd->donatedHydrogen->getGlobalIndex();
859 // no need to register, just look up the index:
860 index = GIDtoH_[frame][hInd];
861
862 dPos = hbd->donorAtom->getPos();
863 hPos = hbd->donatedHydrogen->getPos();
864
865 if (isHBond(dPos, hPos, aPos)) {
866 aInd = hba->getGlobalIndex();
867 registerHydrogenBond(frame, index, hInd, aInd);
868 pos = hPos;
869 RealType shortest = HONKING_LARGE_VALUE;
870
871 // loop over selection 3 to find closest atom in selection 3:
872 for (sd3 = seleMan3_.beginSelected(isd3); sd3 != NULL;
873 sd3 = seleMan3_.nextSelected(isd3)) {
874 vec = pos - sd3->getPos();
875
876 if (usePeriodicBoundaryConditions_)
877 currentSnapshot_->wrapVector(vec);
878
879 r = vec.length();
880
881 if (r < shortest) shortest = r;
882 }
883
884 int whichBin = int(shortest / deltaR_);
885 if (whichBin < int(nRBins_)) { rbin_[frame][index] = whichBin; }
886 }
887 }
888 }
889 }
890 }
891 }
892
893 int HBondJumpR::registerHydrogen(int frame, int hIndex) {
894 int index;
895
896 // If this hydrogen wasn't already registered, register it:
897 if (GIDtoH_[frame][hIndex] == -1) {
898 index = hydrogen_[frame].size();
899
900 GIDtoH_[frame][hIndex] = index;
901 hydrogen_[frame].push_back(hIndex);
902 acceptor_[frame].push_back(-1);
903 selected_[frame].push_back(false);
904 rbin_[frame].push_back(-1);
905
906 if (frame == 0) {
907 lastAcceptor_[frame].push_back(-1);
908 acceptorStartFrame_[frame].push_back(frame);
909 } else {
910 // Copy the last acceptor.
911 int prevIndex = GIDtoH_[frame - 1][hIndex];
912 lastAcceptor_[frame].push_back(lastAcceptor_[frame - 1][prevIndex]);
913 acceptorStartFrame_[frame].push_back(
914 acceptorStartFrame_[frame - 1][prevIndex]);
915 }
916 } else {
917 // This hydrogen was already registered. Just return the index:
918 index = GIDtoH_[frame][hIndex];
919 }
920 return index;
921 }
922
923 void HBondJumpR::correlation() {
924 std::vector<int> s1;
925 std::vector<int>::iterator i1;
926 int index1, index2, gid, aInd1, aInd2;
927
928 for (int i = 0; i < nFrames_; ++i) {
929 RealType time1 = times_[i];
930 s1 = hydrogen_[i];
931
932 for (int j = i; j < nFrames_; ++j) {
933 // Perform a sanity check on the actual configuration times to
934 // make sure the configurations are spaced the same amount the
935 // sample time said they were spaced:
936
937 RealType time2 = times_[j];
938
939 if (fabs((time2 - time1) - (j - i) * deltaTime_) > 1.0e-4) {
940 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
941 "HBondJump::correlation Error: sampleTime (%f)\n"
942 "\tin %s does not match actual time-spacing between\n"
943 "\tconfigurations %d (t = %f) and %d (t = %f).\n",
944 deltaTime_, dumpFilename_.c_str(), i, time1, j, time2);
945 painCave.isFatal = 1;
946 simError();
947 }
948
949 int timeBin = int((time2 - time1) / deltaTime_ + 0.5);
950
951 // loop over the Hydrogens found in frame i:
952
953 for (i1 = s1.begin(); i1 != s1.end(); ++i1) {
954 // gid is the global ID of Hydrogen index1 in frame i:
955 gid = *i1;
956 index1 = GIDtoH_[i][gid];
957
958 // find matching hydrogen in frame j:
959 index2 = GIDtoH_[j][gid];
960
961 if (selected_[i][index1]) {
962 if (int rBin {rbin_[i][index1]}; rBin != -1) {
963 counts_[timeBin][rBin]++;
964
965 if (acceptor_[i][index1] == -1) {
966 aInd1 = lastAcceptor_[i][index1];
967 } else {
968 aInd1 = acceptor_[i][index1];
969 }
970
971 if (acceptor_[j][index2] == -1) {
972 aInd2 = lastAcceptor_[j][index2];
973 } else {
974 aInd2 = acceptor_[j][index2];
975 }
976
977 // aInd1 = acceptor_[i][index1];
978 // aInd2 = acceptor_[j][index2];
979
980 if (aInd1 != aInd2) {
981 // different acceptor so nA(0) . nB(t) = 1
982 histogram_[timeBin][rBin] += 1;
983 } else {
984 // same acceptor, but we need to look at the start frames
985 // for these H-bonds to make sure it is the same H-bond:
986 if (acceptorStartFrame_[i][index1] !=
987 acceptorStartFrame_[j][index2]) {
988 // different start frame, so this is considered a
989 // different H-bond:
990 histogram_[timeBin][rBin] += 1;
991 } else {
992 // same start frame, so this is considered the same H-bond:
993 histogram_[timeBin][rBin] += 0;
994 }
995 }
996 }
997 }
998 }
999 }
1000 }
1001 }
1002
1003 void HBondJumpR::postCorrelate() {
1004 for (unsigned int i = 0; i < nTimeBins_; ++i) {
1005 for (unsigned int j = 0; j < nRBins_; ++j) {
1006 if (counts_[i][j] > 0) {
1007 histogram_[i][j] /= counts_[i][j];
1008 } else {
1009 histogram_[i][j] = 0;
1010 }
1011 histogram_[i][j] = 1.0 - histogram_[i][j];
1012 }
1013 }
1014 }
1015
1016 void HBondJumpR::writeCorrelate() {
1017 ofstream ofs(outputFilename_.c_str());
1018
1019 if (ofs.is_open()) {
1020 Revision r;
1021
1022 ofs << "# " << getCorrFuncType() << "\n";
1023 ofs << "# OpenMD " << r.getFullRevision() << "\n";
1024 ofs << "# " << r.getBuildDate() << "\n";
1025 ofs << "# selection script1: \"" << selectionScript1_;
1026 ofs << "\"\tselection script2: \"" << selectionScript2_ << "\"\n";
1027 if (!paramString_.empty())
1028 ofs << "# parameters: " << paramString_ << "\n";
1029 ofs << "#time\tcorrVal\n";
1030
1031 for (unsigned int i = 0; i < nTimeBins_; ++i) {
1032 ofs << times_[i] - times_[0];
1033
1034 for (unsigned int j = 0; j < nRBins_; ++j) {
1035 ofs << "\t" << histogram_[i][j];
1036 }
1037 ofs << "\n";
1038 }
1039
1040 } else {
1041 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
1042 "HBondJumpR::writeCorrelate Error: fail to open %s\n",
1043 outputFilename_.c_str());
1044 painCave.isFatal = 1;
1045 simError();
1046 }
1047
1048 ofs.close();
1049 }
1050
1051} // namespace OpenMD
OpenMD::Vector::length
Real length()
Returns the length of this vector.
Definition Vector.hpp:393
OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Definition ActionCorrFunc.cpp:60
OpenMD::dot
Real dot(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the dot product of two DynamicVectors.
Definition DynamicVector.hpp:438
OpenMD::getPrefix
std::string getPrefix(const std::string &str)
Definition StringUtils.cpp:169