OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
HBondGeometric.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/HBondGeometric.hpp"
46
47#include <vector>
48
49#include "io/DumpReader.hpp"
51#include "utils/Constants.hpp"
52#include "utils/simError.h"
53
54namespace OpenMD {
55
56 HBondGeometric::HBondGeometric(SimInfo* info, const std::string& filename,
57 const std::string& sele1,
58 const std::string& sele2, double rCut,
59 double thetaCut, int nbins) :
60 StaticAnalyser(info, filename, nbins),
61 selectionScript1_(sele1), seleMan1_(info), evaluator1_(info),
62 selectionScript2_(sele2), seleMan2_(info), evaluator2_(info) {
63 setOutputName(getPrefix(filename) + ".hbg");
64
65 ff_ = info_->getForceField();
66
67 evaluator1_.loadScriptString(sele1);
68 if (!evaluator1_.isDynamic()) {
69 seleMan1_.setSelectionSet(evaluator1_.evaluate());
70 }
71 evaluator2_.loadScriptString(sele2);
72 if (!evaluator2_.isDynamic()) {
73 seleMan2_.setSelectionSet(evaluator2_.evaluate());
74 }
75
76 // Set up cutoff values:
77 nBins_ = nbins;
78 rCut_ = rCut;
79 thetaCut_ = thetaCut;
80
81 nHBonds_.resize(nBins_);
82 nDonor_.resize(nBins_);
83 nAcceptor_.resize(nBins_);
84
85 initializeHistogram();
86 }
87
88 void HBondGeometric::initializeHistogram() {
89 std::fill(nHBonds_.begin(), nHBonds_.end(), 0);
90 std::fill(nDonor_.begin(), nDonor_.end(), 0);
91 std::fill(nAcceptor_.begin(), nAcceptor_.end(), 0);
92 nSelected_ = 0;
93 }
94
95 void HBondGeometric::process() {
96 Molecule* mol1;
97 Molecule* mol2;
98 Molecule::HBondDonor* hbd1;
99 Molecule::HBondDonor* hbd2;
100 std::vector<Molecule::HBondDonor*>::iterator hbdi;
101 std::vector<Molecule::HBondDonor*>::iterator hbdj;
102 std::vector<Atom*>::iterator hbai;
103 std::vector<Atom*>::iterator hbaj;
104 Atom* hba1;
105 Atom* hba2;
106 Vector3d dPos;
107 Vector3d aPos;
108 Vector3d hPos;
109 Vector3d DH;
110 Vector3d DA;
111 RealType DAdist, DHdist, theta, ctheta;
112 int ii, jj;
113 int nHB, nA, nD;
114
115 DumpReader reader(info_, dumpFilename_);
116 int nFrames = reader.getNFrames();
117 frameCounter_ = 0;
118
119 for (int istep = 0; istep < nFrames; istep += step_) {
120 reader.readFrame(istep);
121 frameCounter_++;
122 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
123
124 if (evaluator1_.isDynamic()) {
125 seleMan1_.setSelectionSet(evaluator1_.evaluate());
126 }
127 if (evaluator2_.isDynamic()) {
128 seleMan2_.setSelectionSet(evaluator2_.evaluate());
129 }
130
131 for (mol1 = seleMan1_.beginSelectedMolecule(ii); mol1 != NULL;
132 mol1 = seleMan1_.nextSelectedMolecule(ii)) {
133 // We're collecting statistics on the molecules in selection 1:
134 nHB = 0;
135 nA = 0;
136 nD = 0;
137
138 for (mol2 = seleMan2_.beginSelectedMolecule(jj); mol2 != NULL;
139 mol2 = seleMan2_.nextSelectedMolecule(jj)) {
140 // loop over the possible donors in molecule 1:
141 for (hbd1 = mol1->beginHBondDonor(hbdi); hbd1 != NULL;
142 hbd1 = mol1->nextHBondDonor(hbdi)) {
143 dPos = hbd1->donorAtom->getPos();
144 hPos = hbd1->donatedHydrogen->getPos();
145 DH = hPos - dPos;
146 currentSnapshot_->wrapVector(DH);
147 DHdist = DH.length();
148
149 // loop over the possible acceptors in molecule 2:
150 for (hba2 = mol2->beginHBondAcceptor(hbaj); hba2 != NULL;
151 hba2 = mol2->nextHBondAcceptor(hbaj)) {
152 aPos = hba2->getPos();
153 DA = aPos - dPos;
154 currentSnapshot_->wrapVector(DA);
155 DAdist = DA.length();
156
157 // Distance criteria: are the donor and acceptor atoms
158 // close enough?
159 if (DAdist < rCut_) {
160 ctheta = dot(DH, DA) / (DHdist * DAdist);
161 theta = acos(ctheta) * 180.0 / Constants::PI;
162
163 // Angle criteria: are the D-H and D-A and vectors close?
164 if (theta < thetaCut_) {
165 // molecule 1 is a Hbond donor:
166 nHB++;
167 nD++;
168 }
169 }
170 }
171 }
172
173 // now loop over the possible acceptors in molecule 1:
174 for (hba1 = mol1->beginHBondAcceptor(hbai); hba1 != NULL;
175 hba1 = mol1->nextHBondAcceptor(hbai)) {
176 aPos = hba1->getPos();
177
178 // loop over the possible donors in molecule 2:
179 for (hbd2 = mol2->beginHBondDonor(hbdj); hbd2 != NULL;
180 hbd2 = mol2->nextHBondDonor(hbdj)) {
181 dPos = hbd2->donorAtom->getPos();
182
183 DA = aPos - dPos;
184 currentSnapshot_->wrapVector(DA);
185 DAdist = DA.length();
186
187 // Distance criteria: are the donor and acceptor atoms
188 // close enough?
189 if (DAdist < rCut_) {
190 hPos = hbd2->donatedHydrogen->getPos();
191 DH = hPos - dPos;
192 currentSnapshot_->wrapVector(DH);
193 DHdist = DH.length();
194 ctheta = dot(DH, DA) / (DHdist * DAdist);
195 theta = acos(ctheta) * 180.0 / Constants::PI;
196 // Angle criteria: are the D-H and D-A and vectors close?
197 if (theta < thetaCut_) {
198 // molecule 1 is a Hbond acceptor:
199 nHB++;
200 nA++;
201 }
202 }
203 }
204 }
205 }
206 collectHistogram(nHB, nA, nD);
207 }
208 }
209 writeHistogram();
210 }
211
212 void HBondGeometric::collectHistogram(int nHB, int nA, int nD) {
213 nHBonds_[nHB] += 1;
214 nAcceptor_[nA] += 1;
215 nDonor_[nD] += 1;
216 nSelected_++;
217 }
218
219 void HBondGeometric::writeHistogram() {
220 std::ofstream osq(getOutputFileName().c_str());
221
222 if (osq.is_open()) {
223 osq << "# HydrogenBonding Statistics\n";
224 osq << "# selection1: (" << selectionScript1_ << ")"
225 << "\tselection2: (" << selectionScript2_ << ")\n";
226 osq << "# molecules in selection1: " << nSelected_ << "\n";
227 osq << "# "
228 "nHBonds\tnAcceptor\tnDonor\tp(nHBonds)\tp(nAcceptor)\tp(nDonor)"
229 "\n";
230 // Normalize by number of frames and write it out:
231 for (int i = 0; i < nBins_; ++i) {
232 osq << i;
233 osq << "\t" << nHBonds_[i];
234 osq << "\t" << nAcceptor_[i];
235 osq << "\t" << nDonor_[i];
236 osq << "\t" << (RealType)(nHBonds_[i]) / nSelected_;
237 osq << "\t" << (RealType)(nAcceptor_[i]) / nSelected_;
238 osq << "\t" << (RealType)(nDonor_[i]) / nSelected_;
239 osq << "\n";
240 }
241 osq.close();
242
243 } else {
244 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
245 "HBondGeometric: unable to open %s\n",
246 (getOutputFileName() + "q").c_str());
247 painCave.isFatal = 1;
248 simError();
249 }
250 }
251} // 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)