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