53 |
|
#include "primitives/Molecule.hpp" |
54 |
|
#include "utils/MemoryUtils.hpp" |
55 |
|
#include "utils/simError.h" |
56 |
< |
#include "utils/ElementsTable.hpp" |
56 |
> |
#include "utils/StringUtils.hpp" |
57 |
|
|
58 |
|
namespace OpenMD { |
59 |
|
Molecule::Molecule(int stampId, int globalIndex, const std::string& molName, |
196 |
|
// get the chain of base types for this atom type: |
197 |
|
std::vector<AtomType*> ayb = at->allYourBase(); |
198 |
|
// use the last type in the chain of base types for the name: |
199 |
< |
std::string bn = ayb[ayb.size()-1]->getName(); |
199 |
> |
std::string bn = UpperCase(ayb[ayb.size()-1]->getName()); |
200 |
|
|
201 |
< |
int obanum = etab.GetAtomicNum(bn.c_str()); |
202 |
< |
if (obanum != 0) { |
203 |
< |
RealType eneg = etab.GetElectroNeg(obanum); |
204 |
< |
if (eneg > 3.01) { |
201 |
> |
if (bn.compare("O")==0 || bn.compare("N")==0 |
202 |
> |
|| bn.compare("F")==0) |
203 |
|
hBondAcceptors_.push_back( atom ); |
204 |
< |
} |
207 |
< |
} |
204 |
> |
|
205 |
|
} |
209 |
– |
|
210 |
– |
// find electronegative atoms that are either bonded to hydrogens or are |
211 |
– |
// present in the same rigid bodies: |
206 |
|
|
207 |
+ |
// find electronegative atoms that are either bonded to |
208 |
+ |
// hydrogens or are present in the same rigid bodies: |
209 |
+ |
|
210 |
|
for (bond = beginBond(bi); bond != NULL; bond = nextBond(bi)) { |
211 |
|
Atom* atom1 = bond->getAtomA(); |
212 |
|
Atom* atom2 = bond->getAtomB(); |
216 |
|
std::vector<AtomType*> ayb1 = at1->allYourBase(); |
217 |
|
std::vector<AtomType*> ayb2 = at2->allYourBase(); |
218 |
|
// use the last type in the chain of base types for the name: |
219 |
< |
std::string bn1 = ayb1[ayb1.size()-1]->getName(); |
220 |
< |
std::string bn2 = ayb2[ayb2.size()-1]->getName(); |
221 |
< |
int obanum1 = etab.GetAtomicNum(bn1.c_str()); |
222 |
< |
int obanum2 = etab.GetAtomicNum(bn2.c_str()); |
223 |
< |
|
224 |
< |
if (obanum1 == 1) { |
225 |
< |
if (obanum2 != 0) { |
226 |
< |
RealType eneg = etab.GetElectroNeg(obanum2); |
227 |
< |
if (eneg > 3.01) { |
228 |
< |
HBondDonor* donor = new HBondDonor(); |
232 |
< |
donor->donorAtom = atom2; |
233 |
< |
donor->donatedHydrogen = atom1; |
234 |
< |
hBondDonors_.push_back( donor ); |
235 |
< |
} |
219 |
> |
std::string bn1 = UpperCase(ayb1[ayb1.size()-1]->getName()); |
220 |
> |
std::string bn2 = UpperCase(ayb2[ayb2.size()-1]->getName()); |
221 |
> |
|
222 |
> |
if (bn1.compare("H")==0) { |
223 |
> |
if (bn2.compare("O")==0 || bn2.compare("N")==0 |
224 |
> |
|| bn2.compare("F")==0) { |
225 |
> |
HBondDonor* donor = new HBondDonor(); |
226 |
> |
donor->donorAtom = atom2; |
227 |
> |
donor->donatedHydrogen = atom1; |
228 |
> |
hBondDonors_.push_back( donor ); |
229 |
|
} |
230 |
|
} |
231 |
< |
if (obanum2 == 1) { |
232 |
< |
if (obanum1 != 0) { |
233 |
< |
RealType eneg = etab.GetElectroNeg(obanum1); |
234 |
< |
if (eneg > 3.01) { |
235 |
< |
HBondDonor* donor = new HBondDonor(); |
236 |
< |
donor->donorAtom = atom1; |
237 |
< |
donor->donatedHydrogen = atom2; |
238 |
< |
hBondDonors_.push_back( donor ); |
239 |
< |
} |
247 |
< |
} |
248 |
< |
} |
231 |
> |
if (bn2.compare("H")==0) { |
232 |
> |
if (bn1.compare("O")==0 || bn1.compare("N")==0 |
233 |
> |
|| bn1.compare("F")==0) { |
234 |
> |
HBondDonor* donor = new HBondDonor(); |
235 |
> |
donor->donorAtom = atom1; |
236 |
> |
donor->donatedHydrogen = atom2; |
237 |
> |
hBondDonors_.push_back( donor ); |
238 |
> |
} |
239 |
> |
} |
240 |
|
} |
241 |
< |
|
242 |
< |
for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) { |
243 |
< |
for(atom1 = rb->beginAtom(ai); atom1 != NULL; atom1 = rb->nextAtom(ai)) { |
241 |
> |
|
242 |
> |
for (rb = beginRigidBody(rbIter); rb != NULL; |
243 |
> |
rb = nextRigidBody(rbIter)) { |
244 |
> |
for(atom1 = rb->beginAtom(ai); atom1 != NULL; |
245 |
> |
atom1 = rb->nextAtom(ai)) { |
246 |
|
AtomType* at1 = atom1->getAtomType(); |
247 |
|
// get the chain of base types for this atom type: |
248 |
|
std::vector<AtomType*> ayb1 = at1->allYourBase(); |
249 |
|
// use the last type in the chain of base types for the name: |
250 |
< |
std::string bn1 = ayb1[ayb1.size()-1]->getName(); |
251 |
< |
int obanum1 = etab.GetAtomicNum(bn1.c_str()); |
252 |
< |
if (obanum1 != 0) { |
253 |
< |
RealType eneg = etab.GetElectroNeg(obanum1); |
254 |
< |
if (eneg > 3.01) { |
255 |
< |
for(atom2 = rb->beginAtom(aj); atom2 != NULL; |
256 |
< |
atom2 = rb->nextAtom(aj)) { |
257 |
< |
AtomType* at2 = atom2->getAtomType(); |
258 |
< |
// get the chain of base types for this atom type: |
259 |
< |
std::vector<AtomType*> ayb2 = at2->allYourBase(); |
260 |
< |
// use the last type in the chain of base types for the name: |
261 |
< |
std::string bn2 = ayb2[ayb2.size()-1]->getName(); |
262 |
< |
int obanum2 = etab.GetAtomicNum(bn2.c_str()); |
263 |
< |
if (obanum2 == 1) { |
264 |
< |
HBondDonor* donor = new HBondDonor(); |
265 |
< |
donor->donorAtom = atom1; |
273 |
< |
donor->donatedHydrogen = atom2; |
274 |
< |
hBondDonors_.push_back( donor ); |
275 |
< |
} |
250 |
> |
std::string bn1 = UpperCase(ayb1[ayb1.size()-1]->getName()); |
251 |
> |
|
252 |
> |
if (bn1.compare("O")==0 || bn1.compare("N")==0 |
253 |
> |
|| bn1.compare("F")==0) { |
254 |
> |
for(atom2 = rb->beginAtom(aj); atom2 != NULL; |
255 |
> |
atom2 = rb->nextAtom(aj)) { |
256 |
> |
AtomType* at2 = atom2->getAtomType(); |
257 |
> |
// get the chain of base types for this atom type: |
258 |
> |
std::vector<AtomType*> ayb2 = at2->allYourBase(); |
259 |
> |
// use the last type in the chain of base types for the name: |
260 |
> |
std::string bn2 = UpperCase(ayb2[ayb2.size()-1]->getName()); |
261 |
> |
if (bn2.compare("H")==0) { |
262 |
> |
HBondDonor* donor = new HBondDonor(); |
263 |
> |
donor->donorAtom = atom1; |
264 |
> |
donor->donatedHydrogen = atom2; |
265 |
> |
hBondDonors_.push_back( donor ); |
266 |
|
} |
267 |
|
} |
268 |
|
} |
269 |
|
} |
270 |
< |
} |
270 |
> |
} |
271 |
|
} |
272 |
< |
|
272 |
> |
|
273 |
|
RealType Molecule::getMass() { |
274 |
|
StuntDouble* sd; |
275 |
|
std::vector<StuntDouble*>::iterator i; |