6 |
|
* redistribute this software in source and binary code form, provided |
7 |
|
* that the following conditions are met: |
8 |
|
* |
9 |
< |
* 1. Acknowledgement of the program authors must be made in any |
10 |
< |
* publication of scientific results based in part on use of the |
11 |
< |
* program. An acceptable form of acknowledgement is citation of |
12 |
< |
* the article in which the program was described (Matthew |
13 |
< |
* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher |
14 |
< |
* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented |
15 |
< |
* Parallel Simulation Engine for Molecular Dynamics," |
16 |
< |
* J. Comput. Chem. 26, pp. 252-271 (2005)) |
17 |
< |
* |
18 |
< |
* 2. Redistributions of source code must retain the above copyright |
9 |
> |
* 1. Redistributions of source code must retain the above copyright |
10 |
|
* notice, this list of conditions and the following disclaimer. |
11 |
|
* |
12 |
< |
* 3. Redistributions in binary form must reproduce the above copyright |
12 |
> |
* 2. Redistributions in binary form must reproduce the above copyright |
13 |
|
* notice, this list of conditions and the following disclaimer in the |
14 |
|
* documentation and/or other materials provided with the |
15 |
|
* distribution. |
28 |
|
* arising out of the use of or inability to use software, even if the |
29 |
|
* University of Notre Dame has been advised of the possibility of |
30 |
|
* such damages. |
31 |
+ |
* |
32 |
+ |
* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your |
33 |
+ |
* research, please cite the appropriate papers when you publish your |
34 |
+ |
* work. Good starting points are: |
35 |
+ |
* |
36 |
+ |
* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). |
37 |
+ |
* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). |
38 |
+ |
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). |
39 |
+ |
* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
40 |
+ |
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
41 |
|
*/ |
42 |
|
|
43 |
|
/** |
44 |
|
* @file MoleculeCreator.cpp |
45 |
|
* @author tlin |
46 |
|
* @date 11/04/2004 |
46 |
– |
* @time 13:44am |
47 |
|
* @version 1.0 |
48 |
|
*/ |
49 |
|
|
54 |
|
#include "brains/MoleculeCreator.hpp" |
55 |
|
#include "primitives/GhostBend.hpp" |
56 |
|
#include "primitives/GhostTorsion.hpp" |
57 |
< |
#include "types/DirectionalAtomType.hpp" |
57 |
> |
#include "types/AtomType.hpp" |
58 |
|
#include "types/FixedBondType.hpp" |
59 |
|
#include "utils/simError.h" |
60 |
|
#include "utils/StringUtils.hpp" |
61 |
|
|
62 |
< |
namespace oopse { |
62 |
> |
namespace OpenMD { |
63 |
|
|
64 |
|
Molecule* MoleculeCreator::createMolecule(ForceField* ff, |
65 |
|
MoleculeStamp *molStamp, |
66 |
|
int stampId, int globalIndex, |
67 |
|
LocalIndexManager* localIndexMan) { |
68 |
< |
Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getName()); |
69 |
< |
|
68 |
> |
Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getName(), |
69 |
> |
molStamp->getRegion() ); |
70 |
> |
|
71 |
|
//create atoms |
72 |
|
Atom* atom; |
73 |
|
AtomStamp* currentAtomStamp; |
97 |
|
|
98 |
|
for (int i = 0; i < nBonds; ++i) { |
99 |
|
currentBondStamp = molStamp->getBondStamp(i); |
100 |
< |
bond = createBond(ff, mol, currentBondStamp); |
100 |
> |
bond = createBond(ff, mol, currentBondStamp, localIndexMan); |
101 |
|
mol->addBond(bond); |
102 |
|
} |
103 |
|
|
107 |
|
int nBends = molStamp->getNBends(); |
108 |
|
for (int i = 0; i < nBends; ++i) { |
109 |
|
currentBendStamp = molStamp->getBendStamp(i); |
110 |
< |
bend = createBend(ff, mol, currentBendStamp); |
110 |
> |
bend = createBend(ff, mol, currentBendStamp, localIndexMan); |
111 |
|
mol->addBend(bend); |
112 |
|
} |
113 |
|
|
117 |
|
int nTorsions = molStamp->getNTorsions(); |
118 |
|
for (int i = 0; i < nTorsions; ++i) { |
119 |
|
currentTorsionStamp = molStamp->getTorsionStamp(i); |
120 |
< |
torsion = createTorsion(ff, mol, currentTorsionStamp); |
120 |
> |
torsion = createTorsion(ff, mol, currentTorsionStamp, localIndexMan); |
121 |
|
mol->addTorsion(torsion); |
122 |
|
} |
123 |
|
|
127 |
|
int nInversions = molStamp->getNInversions(); |
128 |
|
for (int i = 0; i < nInversions; ++i) { |
129 |
|
currentInversionStamp = molStamp->getInversionStamp(i); |
130 |
< |
inversion = createInversion(ff, mol, currentInversionStamp); |
130 |
> |
inversion = createInversion(ff, mol, currentInversionStamp, |
131 |
> |
localIndexMan); |
132 |
|
if (inversion != NULL ) { |
133 |
|
mol->addInversion(inversion); |
134 |
|
} |
140 |
|
int nCutoffGroups = molStamp->getNCutoffGroups(); |
141 |
|
for (int i = 0; i < nCutoffGroups; ++i) { |
142 |
|
currentCutoffGroupStamp = molStamp->getCutoffGroupStamp(i); |
143 |
< |
cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp); |
143 |
> |
cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp, |
144 |
> |
localIndexMan); |
145 |
|
mol->addCutoffGroup(cutoffGroup); |
146 |
|
} |
147 |
|
|
172 |
|
// every single free atom |
173 |
|
|
174 |
|
for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) { |
175 |
< |
cutoffGroup = createCutoffGroup(mol, *fai); |
175 |
> |
cutoffGroup = createCutoffGroup(mol, *fai, localIndexMan); |
176 |
|
mol->addCutoffGroup(cutoffGroup); |
177 |
|
} |
178 |
|
//create constraints |
179 |
|
createConstraintPair(mol); |
180 |
|
createConstraintElem(mol); |
181 |
|
|
182 |
+ |
// Does this molecule stamp define a total constrained charge value? |
183 |
+ |
// If so, let the created molecule know about it. |
184 |
+ |
|
185 |
+ |
if (molStamp->haveConstrainTotalCharge() ) { |
186 |
+ |
mol->setConstrainTotalCharge( molStamp->getConstrainTotalCharge() ); |
187 |
+ |
} |
188 |
+ |
|
189 |
|
//the construction of this molecule is finished |
190 |
|
mol->complete(); |
191 |
|
|
208 |
|
painCave.isFatal = 1; |
209 |
|
simError(); |
210 |
|
} |
211 |
< |
|
211 |
> |
|
212 |
|
//below code still have some kind of hard-coding smell |
213 |
|
if (atomType->isDirectional()){ |
204 |
– |
|
205 |
– |
DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType); |
206 |
– |
|
207 |
– |
if (dAtomType == NULL) { |
208 |
– |
sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType"); |
214 |
|
|
210 |
– |
painCave.isFatal = 1; |
211 |
– |
simError(); |
212 |
– |
} |
213 |
– |
|
215 |
|
DirectionalAtom* dAtom; |
216 |
< |
dAtom = new DirectionalAtom(dAtomType); |
216 |
> |
dAtom = new DirectionalAtom(atomType); |
217 |
|
atom = dAtom; |
218 |
|
} |
219 |
|
else{ |
228 |
|
RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, |
229 |
|
Molecule* mol, |
230 |
|
RigidBodyStamp* rbStamp, |
231 |
< |
LocalIndexManager* localIndexMan) { |
231 |
> |
LocalIndexManager* localIndexMan){ |
232 |
|
Atom* atom; |
233 |
|
int nAtoms; |
234 |
|
Vector3d refCoor; |
252 |
|
//set the local index of this rigid body, global index will be set later |
253 |
|
rb->setLocalIndex(localIndexMan->getNextRigidBodyIndex()); |
254 |
|
|
255 |
< |
//the rule for naming rigidbody MoleculeName_RB_Integer |
256 |
< |
//The first part is the name of the molecule |
257 |
< |
//The second part is alway fixed as "RB" |
258 |
< |
//The third part is the index of the rigidbody defined in meta-data file |
259 |
< |
//For example, Butane_RB_0 is a valid rigid body name of butane molecule |
259 |
< |
/**@todo replace itoa by lexi_cast */ |
260 |
< |
std::string s = OOPSE_itoa(mol->getNRigidBodies(), 10); |
261 |
< |
rb->setType(mol->getType() + "_RB_" + s.c_str()); |
255 |
> |
// The rule for naming a rigidbody is: MoleculeName_RB_Integer |
256 |
> |
// The first part is the name of the molecule |
257 |
> |
// The second part is always fixed as "RB" |
258 |
> |
// The third part is the index of the rigidbody defined in meta-data file |
259 |
> |
// For example, Butane_RB_0 is a valid rigid body name of butane molecule |
260 |
|
|
261 |
+ |
std::string s = OpenMD_itoa(mol->getNRigidBodies(), 10); |
262 |
+ |
rb->setType(mol->getType() + "_RB_" + s.c_str()); |
263 |
|
return rb; |
264 |
|
} |
265 |
|
|
266 |
|
Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol, |
267 |
< |
BondStamp* stamp) { |
267 |
> |
BondStamp* stamp, |
268 |
> |
LocalIndexManager* localIndexMan) { |
269 |
|
BondType* bondType; |
270 |
|
Atom* atomA; |
271 |
|
Atom* atomB; |
285 |
|
painCave.isFatal = 1; |
286 |
|
simError(); |
287 |
|
} |
288 |
< |
return new Bond(atomA, atomB, bondType); |
288 |
> |
Bond* bond = new Bond(atomA, atomB, bondType); |
289 |
> |
|
290 |
> |
//set the local index of this bond, the global index will be set later |
291 |
> |
bond->setLocalIndex(localIndexMan->getNextBondIndex()); |
292 |
> |
|
293 |
> |
// The rule for naming a bond is: MoleculeName_Bond_Integer |
294 |
> |
// The first part is the name of the molecule |
295 |
> |
// The second part is always fixed as "Bond" |
296 |
> |
// The third part is the index of the bond defined in meta-data file |
297 |
> |
// For example, Butane_bond_0 is a valid Bond name in a butane molecule |
298 |
> |
|
299 |
> |
std::string s = OpenMD_itoa(mol->getNBonds(), 10); |
300 |
> |
bond->setName(mol->getType() + "_Bond_" + s.c_str()); |
301 |
> |
return bond; |
302 |
|
} |
303 |
|
|
304 |
|
Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, |
305 |
< |
BendStamp* stamp) { |
305 |
> |
BendStamp* stamp, |
306 |
> |
LocalIndexManager* localIndexMan) { |
307 |
|
Bend* bend = NULL; |
308 |
|
std::vector<int> bendAtoms = stamp->getMembers(); |
309 |
|
if (bendAtoms.size() == 3) { |
318 |
|
atomC->getType().c_str()); |
319 |
|
|
320 |
|
if (bendType == NULL) { |
321 |
< |
sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]", |
321 |
> |
sprintf(painCave.errMsg, |
322 |
> |
"Can not find Matching Bend Type for[%s, %s, %s]", |
323 |
|
atomA->getType().c_str(), |
324 |
|
atomB->getType().c_str(), |
325 |
|
atomC->getType().c_str()); |
343 |
|
BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST"); |
344 |
|
|
345 |
|
if (bendType == NULL) { |
346 |
< |
sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]", |
346 |
> |
sprintf(painCave.errMsg, |
347 |
> |
"Can not find Matching Bend Type for[%s, %s, %s]", |
348 |
|
normalAtom->getType().c_str(), |
349 |
|
ghostAtom->getType().c_str(), |
350 |
|
"GHOST"); |
356 |
|
bend = new GhostBend(normalAtom, ghostAtom, bendType); |
357 |
|
|
358 |
|
} |
359 |
< |
|
359 |
> |
|
360 |
> |
//set the local index of this bend, the global index will be set later |
361 |
> |
bend->setLocalIndex(localIndexMan->getNextBendIndex()); |
362 |
> |
|
363 |
> |
// The rule for naming a bend is: MoleculeName_Bend_Integer |
364 |
> |
// The first part is the name of the molecule |
365 |
> |
// The second part is always fixed as "Bend" |
366 |
> |
// The third part is the index of the bend defined in meta-data file |
367 |
> |
// For example, Butane_Bend_0 is a valid Bend name in a butane molecule |
368 |
> |
|
369 |
> |
std::string s = OpenMD_itoa(mol->getNBends(), 10); |
370 |
> |
bend->setName(mol->getType() + "_Bend_" + s.c_str()); |
371 |
|
return bend; |
372 |
|
} |
373 |
|
|
374 |
|
Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, |
375 |
< |
TorsionStamp* stamp) { |
375 |
> |
TorsionStamp* stamp, |
376 |
> |
LocalIndexManager* localIndexMan) { |
377 |
|
|
378 |
|
Torsion* torsion = NULL; |
379 |
|
std::vector<int> torsionAtoms = stamp->getMembers(); |
394 |
|
atomB->getType(), |
395 |
|
atomC->getType(), |
396 |
|
atomD->getType()); |
368 |
– |
|
397 |
|
if (torsionType == NULL) { |
398 |
< |
sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]", |
398 |
> |
sprintf(painCave.errMsg, |
399 |
> |
"Can not find Matching Torsion Type for[%s, %s, %s, %s]", |
400 |
|
atomA->getType().c_str(), |
401 |
|
atomB->getType().c_str(), |
402 |
|
atomC->getType().c_str(), |
433 |
|
|
434 |
|
torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType); |
435 |
|
} |
436 |
+ |
|
437 |
+ |
//set the local index of this torsion, the global index will be set later |
438 |
+ |
torsion->setLocalIndex(localIndexMan->getNextTorsionIndex()); |
439 |
|
|
440 |
+ |
// The rule for naming a torsion is: MoleculeName_Torsion_Integer |
441 |
+ |
// The first part is the name of the molecule |
442 |
+ |
// The second part is always fixed as "Torsion" |
443 |
+ |
// The third part is the index of the torsion defined in meta-data file |
444 |
+ |
// For example, Butane_Torsion_0 is a valid Torsion name in a |
445 |
+ |
// butane molecule |
446 |
+ |
|
447 |
+ |
std::string s = OpenMD_itoa(mol->getNTorsions(), 10); |
448 |
+ |
torsion->setName(mol->getType() + "_Torsion_" + s.c_str()); |
449 |
|
return torsion; |
450 |
|
} |
451 |
|
|
452 |
|
Inversion* MoleculeCreator::createInversion(ForceField* ff, Molecule* mol, |
453 |
< |
InversionStamp* stamp) { |
453 |
> |
InversionStamp* stamp, |
454 |
> |
LocalIndexManager* localIndexMan) { |
455 |
|
|
456 |
|
Inversion* inversion = NULL; |
457 |
|
int center = stamp->getCenter(); |
481 |
|
atomD->getType().c_str()); |
482 |
|
|
483 |
|
painCave.isFatal = 0; |
484 |
< |
painCave.severity = OOPSE_INFO; |
484 |
> |
painCave.severity = OPENMD_INFO; |
485 |
|
simError(); |
486 |
|
return NULL; |
487 |
|
} else { |
488 |
|
|
489 |
|
inversion = new Inversion(atomA, atomB, atomC, atomD, inversionType); |
490 |
+ |
|
491 |
+ |
// set the local index of this inversion, the global index will |
492 |
+ |
// be set later |
493 |
+ |
inversion->setLocalIndex(localIndexMan->getNextInversionIndex()); |
494 |
+ |
|
495 |
+ |
// The rule for naming an inversion is: MoleculeName_Inversion_Integer |
496 |
+ |
// The first part is the name of the molecule |
497 |
+ |
// The second part is always fixed as "Inversion" |
498 |
+ |
// The third part is the index of the inversion defined in meta-data file |
499 |
+ |
// For example, Benzene_Inversion_0 is a valid Inversion name in a |
500 |
+ |
// Benzene molecule |
501 |
+ |
|
502 |
+ |
std::string s = OpenMD_itoa(mol->getNInversions(), 10); |
503 |
+ |
inversion->setName(mol->getType() + "_Inversion_" + s.c_str()); |
504 |
|
return inversion; |
505 |
|
} |
506 |
|
} |
507 |
|
|
508 |
|
|
509 |
< |
CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) { |
509 |
> |
CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, |
510 |
> |
CutoffGroupStamp* stamp, |
511 |
> |
LocalIndexManager* localIndexMan) { |
512 |
|
int nAtoms; |
513 |
|
CutoffGroup* cg; |
514 |
|
Atom* atom; |
520 |
|
assert(atom); |
521 |
|
cg->addAtom(atom); |
522 |
|
} |
523 |
< |
|
523 |
> |
|
524 |
> |
//set the local index of this cutoffGroup, global index will be set later |
525 |
> |
cg->setLocalIndex(localIndexMan->getNextCutoffGroupIndex()); |
526 |
> |
|
527 |
|
return cg; |
528 |
|
} |
529 |
< |
|
530 |
< |
CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) { |
529 |
> |
|
530 |
> |
CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom, |
531 |
> |
LocalIndexManager* localIndexMan) { |
532 |
|
CutoffGroup* cg; |
533 |
|
cg = new CutoffGroup(); |
534 |
|
cg->addAtom(atom); |
535 |
+ |
|
536 |
+ |
//set the local index of this cutoffGroup, global index will be set later |
537 |
+ |
cg->setLocalIndex(localIndexMan->getNextCutoffGroupIndex()); |
538 |
+ |
|
539 |
|
return cg; |
540 |
|
} |
541 |
|
|
548 |
|
|
549 |
|
BondType* bt = bond->getBondType(); |
550 |
|
|
485 |
– |
//class Parent1 {}; |
486 |
– |
//class Child1 : public Parent {}; |
487 |
– |
//class Child2 : public Parent {}; |
488 |
– |
//Child1* ch1 = new Child1(); |
489 |
– |
//Child2* ch2 = dynamic_cast<Child2*>(ch1); |
490 |
– |
//the dynamic_cast is succeed in above line. A compiler bug? |
491 |
– |
|
551 |
|
if (typeid(FixedBondType) == typeid(*bt)) { |
552 |
|
FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt); |
553 |
|
|