| 56 |
|
const std::string& sele2, |
| 57 |
|
double rCut, double thetaCut, int nbins) : |
| 58 |
|
StaticAnalyser(info, filename), |
| 59 |
< |
selectionScript1_(sele1), evaluator1_(info), seleMan1_(info), |
| 60 |
< |
selectionScript2_(sele2), evaluator2_(info), seleMan2_(info){ |
| 59 |
> |
selectionScript1_(sele1), seleMan1_(info), evaluator1_(info), |
| 60 |
> |
selectionScript2_(sele2), seleMan2_(info), evaluator2_(info) { |
| 61 |
|
|
| 62 |
|
setOutputName(getPrefix(filename) + ".hbg"); |
| 63 |
|
|
| 81 |
|
nHBonds_.resize(nBins_); |
| 82 |
|
nDonor_.resize(nBins_); |
| 83 |
|
nAcceptor_.resize(nBins_); |
| 84 |
+ |
|
| 85 |
+ |
initializeHistogram(); |
| 86 |
|
} |
| 87 |
|
|
| 88 |
|
HBondGeometric::~HBondGeometric() { |
| 98 |
|
nSelected_ = 0; |
| 99 |
|
} |
| 100 |
|
|
| 99 |
– |
|
| 100 |
– |
|
| 101 |
|
void HBondGeometric::process() { |
| 102 |
< |
Molecule* mol; |
| 103 |
< |
StuntDouble* sd1; |
| 104 |
< |
StuntDouble* sd2; |
| 102 |
> |
Molecule* mol1; |
| 103 |
> |
Molecule* mol2; |
| 104 |
|
RigidBody* rb1; |
| 105 |
< |
RigidBody* rb2; |
| 105 |
> |
Molecule::HBondDonor* hbd1; |
| 106 |
> |
Molecule::HBondDonor* hbd2; |
| 107 |
> |
std::vector<Molecule::HBondDonor*>::iterator hbdi; |
| 108 |
> |
std::vector<Molecule::HBondDonor*>::iterator hbdj; |
| 109 |
> |
std::vector<Atom*>::iterator hbai; |
| 110 |
> |
std::vector<Atom*>::iterator hbaj; |
| 111 |
> |
Atom* hba1; |
| 112 |
> |
Atom* hba2; |
| 113 |
|
SimInfo::MoleculeIterator mi; |
| 114 |
|
Molecule::RigidBodyIterator rbIter; |
| 115 |
< |
Molecule::IntegrableObjectIterator ioi; |
| 115 |
> |
Vector3d dPos; |
| 116 |
> |
Vector3d aPos; |
| 117 |
> |
Vector3d hPos; |
| 118 |
> |
Vector3d DH; |
| 119 |
> |
Vector3d DA; |
| 120 |
> |
RealType DAdist, DHdist, theta, ctheta; |
| 121 |
|
int ii, jj; |
| 111 |
– |
std::string rbName; |
| 112 |
– |
std::vector<Atom *> atoms1; |
| 113 |
– |
std::vector<Atom *> atoms2; |
| 114 |
– |
std::vector<Atom *>::iterator ai1; |
| 115 |
– |
std::vector<Atom *>::iterator ai2; |
| 116 |
– |
Vector3d O1pos, O2pos; |
| 117 |
– |
Vector3d H1apos, H1bpos, H2apos, H2bpos; |
| 122 |
|
int nHB, nA, nD; |
| 123 |
|
|
| 124 |
|
DumpReader reader(info_, dumpFilename_); |
| 132 |
|
|
| 133 |
|
// update the positions of atoms which belong to the rigidbodies |
| 134 |
|
|
| 135 |
< |
for (mol = info_->beginMolecule(mi); mol != NULL; |
| 136 |
< |
mol = info_->nextMolecule(mi)) { |
| 137 |
< |
for (rb1 = mol->beginRigidBody(rbIter); rb1 != NULL; |
| 138 |
< |
rb1 = mol->nextRigidBody(rbIter)) { |
| 135 |
> |
for (mol1 = info_->beginMolecule(mi); mol1 != NULL; |
| 136 |
> |
mol1 = info_->nextMolecule(mi)) { |
| 137 |
> |
for (rb1 = mol1->beginRigidBody(rbIter); rb1 != NULL; |
| 138 |
> |
rb1 = mol1->nextRigidBody(rbIter)) { |
| 139 |
|
rb1->updateAtoms(); |
| 140 |
|
} |
| 141 |
|
} |
| 147 |
|
seleMan2_.setSelectionSet(evaluator2_.evaluate()); |
| 148 |
|
} |
| 149 |
|
|
| 150 |
< |
for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL; sd1 = seleMan1_.nextSelected(ii)) { |
| 151 |
< |
if (sd1->isRigidBody()) { |
| 148 |
< |
rb1 = dynamic_cast<RigidBody*>(sd1); |
| 149 |
< |
atoms1 = rb1->getAtoms(); |
| 150 |
< |
|
| 151 |
< |
int nH = 0; |
| 152 |
< |
int nO = 0; |
| 153 |
< |
|
| 154 |
< |
for (ai1 = atoms1.begin(); ai1 != atoms1.end(); ++ai1) { |
| 155 |
< |
std::string atName = (*ai1)->getType(); |
| 156 |
< |
// query the force field for the AtomType associated with this |
| 157 |
< |
// atomTypeName: |
| 158 |
< |
AtomType* at = ff_->getAtomType(atName); |
| 159 |
< |
// get the chain of base types for this atom type: |
| 160 |
< |
std::vector<AtomType*> ayb = at->allYourBase(); |
| 161 |
< |
// use the last type in the chain of base types for the name: |
| 162 |
< |
std::string bn = ayb[ayb.size()-1]->getName(); |
| 163 |
< |
|
| 164 |
< |
bool isH = bn.compare("H") == 0 ? true : false; |
| 165 |
< |
bool isO = bn.compare("O") == 0 ? true : false; |
| 166 |
< |
|
| 167 |
< |
if (isO && nO == 0) { |
| 168 |
< |
O1pos = (*ai1)->getPos(); |
| 169 |
< |
nO++; |
| 170 |
< |
} |
| 171 |
< |
if (isH) { |
| 172 |
< |
if (nH == 0) { |
| 173 |
< |
H1apos = (*ai1)->getPos(); |
| 174 |
< |
} |
| 175 |
< |
if (nH == 1) { |
| 176 |
< |
H1bpos = (*ai1)->getPos(); |
| 177 |
< |
} |
| 178 |
< |
nH++; |
| 179 |
< |
} |
| 180 |
< |
} |
| 181 |
< |
} |
| 150 |
> |
for (mol1 = seleMan1_.beginSelectedMolecule(ii); |
| 151 |
> |
mol1 != NULL; mol1 = seleMan1_.nextSelectedMolecule(ii)) { |
| 152 |
|
|
| 153 |
< |
|
| 153 |
> |
// We're collecting statistics on the molecules in selection 1: |
| 154 |
|
nHB = 0; |
| 155 |
|
nA = 0; |
| 156 |
|
nD = 0; |
| 157 |
|
|
| 158 |
< |
for (sd2 = seleMan2_.beginSelected(jj); sd2 != NULL; sd2 = seleMan2_.nextSelected(jj)) { |
| 158 |
> |
for (mol2 = seleMan2_.beginSelectedMolecule(jj); |
| 159 |
> |
mol2 != NULL; mol2 = seleMan2_.nextSelectedMolecule(jj)) { |
| 160 |
> |
|
| 161 |
> |
// loop over the possible donors in molecule 1: |
| 162 |
> |
for (hbd1 = mol1->beginHBondDonor(hbdi); hbd1 != NULL; |
| 163 |
> |
hbd1 = mol1->nextHBondDonor(hbdi)) { |
| 164 |
> |
dPos = hbd1->donorAtom->getPos(); |
| 165 |
> |
hPos = hbd1->donatedHydrogen->getPos(); |
| 166 |
> |
DH = hPos - dPos; |
| 167 |
> |
currentSnapshot_->wrapVector(DH); |
| 168 |
> |
DHdist = DH.length(); |
| 169 |
|
|
| 170 |
< |
if (sd1 == sd2) continue; |
| 171 |
< |
|
| 172 |
< |
if (sd2->isRigidBody()) { |
| 173 |
< |
rb2 = dynamic_cast<RigidBody*>(sd2); |
| 174 |
< |
atoms2 = rb2->getAtoms(); |
| 175 |
< |
|
| 176 |
< |
int nH = 0; |
| 197 |
< |
int nO = 0; |
| 198 |
< |
|
| 199 |
< |
for (ai2 = atoms2.begin(); ai2 != atoms2.end(); ++ai2) { |
| 200 |
< |
std::string atName = (*ai2)->getType(); |
| 201 |
< |
// query the force field for the AtomType associated with this |
| 202 |
< |
// atomTypeName: |
| 203 |
< |
AtomType* at = ff_->getAtomType(atName); |
| 204 |
< |
// get the chain of base types for this atom type: |
| 205 |
< |
std::vector<AtomType*> ayb = at->allYourBase(); |
| 206 |
< |
// use the last type in the chain of base types for the name: |
| 207 |
< |
std::string bn = ayb[ayb.size()-1]->getName(); |
| 208 |
< |
|
| 209 |
< |
bool isH = bn.compare("H") == 0 ? true : false; |
| 210 |
< |
bool isO = bn.compare("O") == 0 ? true : false; |
| 170 |
> |
// loop over the possible acceptors in molecule 2: |
| 171 |
> |
for (hba2 = mol2->beginHBondAcceptor(hbaj); hba2 != NULL; |
| 172 |
> |
hba2 = mol2->nextHBondAcceptor(hbaj)) { |
| 173 |
> |
aPos = hba2->getPos(); |
| 174 |
> |
DA = aPos - dPos; |
| 175 |
> |
currentSnapshot_->wrapVector(DA); |
| 176 |
> |
DAdist = DA.length(); |
| 177 |
|
|
| 178 |
< |
if (isO && nO == 0) { |
| 179 |
< |
O2pos = (*ai2)->getPos(); |
| 180 |
< |
nO++; |
| 181 |
< |
} |
| 182 |
< |
if (isH) { |
| 183 |
< |
if (nH == 0) { |
| 184 |
< |
H2apos = (*ai2)->getPos(); |
| 178 |
> |
// Distance criteria: are the donor and acceptor atoms |
| 179 |
> |
// close enough? |
| 180 |
> |
if (DAdist < rCut_) { |
| 181 |
> |
|
| 182 |
> |
ctheta = dot(DH, DA) / (DHdist * DAdist); |
| 183 |
> |
theta = acos(ctheta) * 180.0 / M_PI; |
| 184 |
> |
|
| 185 |
> |
// Angle criteria: are the D-H and D-A and vectors close? |
| 186 |
> |
if (theta < thetaCut_) { |
| 187 |
> |
// molecule 1 is a Hbond donor: |
| 188 |
> |
nHB++; |
| 189 |
> |
nD++; |
| 190 |
|
} |
| 191 |
< |
if (nH == 1) { |
| 192 |
< |
H2bpos = (*ai2)->getPos(); |
| 193 |
< |
} |
| 194 |
< |
nH++; |
| 195 |
< |
} |
| 196 |
< |
} |
| 191 |
> |
} |
| 192 |
> |
} |
| 193 |
> |
} |
| 194 |
> |
|
| 195 |
> |
// now loop over the possible acceptors in molecule 1: |
| 196 |
> |
for (hba1 = mol1->beginHBondAcceptor(hbai); hba1 != NULL; |
| 197 |
> |
hba1 = mol1->nextHBondAcceptor(hbai)) { |
| 198 |
> |
aPos = hba1->getPos(); |
| 199 |
|
|
| 200 |
< |
// Do our testing: |
| 201 |
< |
Vector3d Odiff = O2pos - O1pos; |
| 202 |
< |
currentSnapshot_->wrapVector(Odiff); |
| 203 |
< |
RealType Odist = Odiff.length(); |
| 231 |
< |
if (Odist < rCut_) { |
| 232 |
< |
// OH vectors: |
| 233 |
< |
Vector3d HO1a = H1apos - O1pos; |
| 234 |
< |
Vector3d HO1b = H1bpos - O1pos; |
| 235 |
< |
Vector3d HO2a = H2apos - O2pos; |
| 236 |
< |
Vector3d HO2b = H2bpos - O2pos; |
| 237 |
< |
// wrapped in case a molecule is split across boundaries: |
| 238 |
< |
currentSnapshot_->wrapVector(HO1a); |
| 239 |
< |
currentSnapshot_->wrapVector(HO1b); |
| 240 |
< |
currentSnapshot_->wrapVector(HO2a); |
| 241 |
< |
currentSnapshot_->wrapVector(HO2a); |
| 242 |
< |
// cos thetas: |
| 243 |
< |
RealType ctheta1a = dot(HO1a, Odiff) / (Odist * HO1a.length()); |
| 244 |
< |
RealType ctheta1b = dot(HO1b, Odiff) / (Odist * HO1b.length()); |
| 245 |
< |
RealType ctheta2a = dot(HO2a, -Odiff) / (Odist * HO2a.length()); |
| 246 |
< |
RealType ctheta2b = dot(HO2b, -Odiff) / (Odist * HO2b.length()); |
| 200 |
> |
// loop over the possible donors in molecule 2: |
| 201 |
> |
for (hbd2 = mol2->beginHBondDonor(hbdj); hbd2 != NULL; |
| 202 |
> |
hbd2 = mol2->nextHBondDonor(hbdj)) { |
| 203 |
> |
dPos = hbd2->donorAtom->getPos(); |
| 204 |
|
|
| 205 |
< |
RealType theta1a = acos(ctheta1a) * 180.0 / M_PI; |
| 206 |
< |
RealType theta1b = acos(ctheta1b) * 180.0 / M_PI; |
| 207 |
< |
RealType theta2a = acos(ctheta2a) * 180.0 / M_PI; |
| 208 |
< |
RealType theta2b = acos(ctheta2b) * 180.0 / M_PI; |
| 209 |
< |
|
| 210 |
< |
if (theta1a < thetaCut_) { |
| 211 |
< |
// molecule 1 is a Hbond donor: |
| 212 |
< |
nHB++; |
| 213 |
< |
nD++; |
| 205 |
> |
DA = aPos - dPos; |
| 206 |
> |
currentSnapshot_->wrapVector(DA); |
| 207 |
> |
DAdist = DA.length(); |
| 208 |
> |
|
| 209 |
> |
// Distance criteria: are the donor and acceptor atoms |
| 210 |
> |
// close enough? |
| 211 |
> |
if (DAdist < rCut_) { |
| 212 |
> |
hPos = hbd2->donatedHydrogen->getPos(); |
| 213 |
> |
DH = hPos - dPos; |
| 214 |
> |
currentSnapshot_->wrapVector(DH); |
| 215 |
> |
DHdist = DH.length(); |
| 216 |
> |
ctheta = dot(DH, DA) / (DHdist * DAdist); |
| 217 |
> |
theta = acos(ctheta) * 180.0 / M_PI; |
| 218 |
> |
// Angle criteria: are the D-H and D-A and vectors close? |
| 219 |
> |
if (theta < thetaCut_) { |
| 220 |
> |
// molecule 1 is a Hbond acceptor: |
| 221 |
> |
nHB++; |
| 222 |
> |
nA++; |
| 223 |
> |
} |
| 224 |
|
} |
| 225 |
< |
if (theta1b < thetaCut_) { |
| 259 |
< |
// molecule 1 is a Hbond donor: |
| 260 |
< |
nHB++; |
| 261 |
< |
nD++; |
| 262 |
< |
} |
| 263 |
< |
if (theta2a < thetaCut_) { |
| 264 |
< |
// molecule 1 is a Hbond acceptor: |
| 265 |
< |
nHB++; |
| 266 |
< |
nA++; |
| 267 |
< |
} |
| 268 |
< |
if (theta2b < thetaCut_) { |
| 269 |
< |
// molecule 1 is a Hbond acceptor: |
| 270 |
< |
nHB++; |
| 271 |
< |
nA++; |
| 272 |
< |
} |
| 273 |
< |
} |
| 225 |
> |
} |
| 226 |
|
} |
| 227 |
< |
} |
| 227 |
> |
} |
| 228 |
|
collectHistogram(nHB, nA, nD); |
| 229 |
|
} |
| 230 |
|
} |
| 243 |
|
void HBondGeometric::writeHistogram() { |
| 244 |
|
|
| 245 |
|
std::ofstream osq(getOutputFileName().c_str()); |
| 294 |
– |
cerr << "nSelected = " << nSelected_ << "\n"; |
| 246 |
|
|
| 247 |
|
if (osq.is_open()) { |
| 248 |
|
|
| 249 |
|
osq << "# HydrogenBonding Statistics\n"; |
| 250 |
|
osq << "# selection1: (" << selectionScript1_ << ")" |
| 251 |
|
<< "\tselection2: (" << selectionScript2_ << ")\n"; |
| 252 |
< |
osq << "# p(nHBonds)\tp(nAcceptor)\tp(nDonor)\n"; |
| 252 |
> |
osq << "# molecules in selection1: " << nSelected_ << "\n"; |
| 253 |
> |
osq << "# nHBonds\tnAcceptor\tnDonor\tp(nHBonds)\tp(nAcceptor)\tp(nDonor)\n"; |
| 254 |
|
// Normalize by number of frames and write it out: |
| 255 |
|
for (int i = 0; i < nBins_; ++i) { |
| 256 |
|
osq << i; |
| 257 |
+ |
osq << "\t" << nHBonds_[i]; |
| 258 |
+ |
osq << "\t" << nAcceptor_[i]; |
| 259 |
+ |
osq << "\t" << nDonor_[i]; |
| 260 |
|
osq << "\t" << (RealType) (nHBonds_[i]) / nSelected_; |
| 261 |
|
osq << "\t" << (RealType) (nAcceptor_[i]) / nSelected_; |
| 262 |
|
osq << "\t" << (RealType) (nDonor_[i]) / nSelected_; |
