ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/brains/MoleculeCreator.cpp
Revision: 1540
Committed: Mon Jan 17 21:34:36 2011 UTC (14 years, 3 months ago) by gezelter
File size: 18115 byte(s)
Log Message:
changes for new parallel architecture

File Contents

# Content
1 /*
2 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 *
4 * The University of Notre Dame grants you ("Licensee") a
5 * non-exclusive, royalty free, license to use, modify and
6 * redistribute this software in source and binary code form, provided
7 * that the following conditions are met:
8 *
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 *
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.
16 *
17 * This software is provided "AS IS," without a warranty of any
18 * kind. All express or implied conditions, representations and
19 * warranties, including any implied warranty of merchantability,
20 * fitness for a particular purpose or non-infringement, are hereby
21 * excluded. The University of Notre Dame and its licensors shall not
22 * be liable for any damages suffered by licensee as a result of
23 * using, modifying or distributing the software or its
24 * derivatives. In no event will the University of Notre Dame or its
25 * licensors be liable for any lost revenue, profit or data, or for
26 * direct, indirect, special, consequential, incidental or punitive
27 * damages, however caused and regardless of the theory of liability,
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, 24107 (2008).
39 * [4] Vardeman & Gezelter, in progress (2009).
40 */
41
42 /**
43 * @file MoleculeCreator.cpp
44 * @author tlin
45 * @date 11/04/2004
46 * @time 13:44am
47 * @version 1.0
48 */
49
50 #include <cassert>
51 #include <typeinfo>
52 #include <set>
53
54 #include "brains/MoleculeCreator.hpp"
55 #include "primitives/GhostBend.hpp"
56 #include "primitives/GhostTorsion.hpp"
57 #include "types/DirectionalAtomType.hpp"
58 #include "types/FixedBondType.hpp"
59 #include "utils/simError.h"
60 #include "utils/StringUtils.hpp"
61
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
70 //create atoms
71 Atom* atom;
72 AtomStamp* currentAtomStamp;
73 int nAtom = molStamp->getNAtoms();
74 for (int i = 0; i < nAtom; ++i) {
75 currentAtomStamp = molStamp->getAtomStamp(i);
76 atom = createAtom(ff, mol, currentAtomStamp, localIndexMan);
77 mol->addAtom(atom);
78 }
79
80 //create rigidbodies
81 RigidBody* rb;
82 RigidBodyStamp * currentRigidBodyStamp;
83 int nRigidbodies = molStamp->getNRigidBodies();
84
85 for (int i = 0; i < nRigidbodies; ++i) {
86 currentRigidBodyStamp = molStamp->getRigidBodyStamp(i);
87 rb = createRigidBody(molStamp, mol, currentRigidBodyStamp,
88 localIndexMan);
89 mol->addRigidBody(rb);
90 }
91
92 //create bonds
93 Bond* bond;
94 BondStamp* currentBondStamp;
95 int nBonds = molStamp->getNBonds();
96
97 for (int i = 0; i < nBonds; ++i) {
98 currentBondStamp = molStamp->getBondStamp(i);
99 bond = createBond(ff, mol, currentBondStamp);
100 mol->addBond(bond);
101 }
102
103 //create bends
104 Bend* bend;
105 BendStamp* currentBendStamp;
106 int nBends = molStamp->getNBends();
107 for (int i = 0; i < nBends; ++i) {
108 currentBendStamp = molStamp->getBendStamp(i);
109 bend = createBend(ff, mol, currentBendStamp);
110 mol->addBend(bend);
111 }
112
113 //create torsions
114 Torsion* torsion;
115 TorsionStamp* currentTorsionStamp;
116 int nTorsions = molStamp->getNTorsions();
117 for (int i = 0; i < nTorsions; ++i) {
118 currentTorsionStamp = molStamp->getTorsionStamp(i);
119 torsion = createTorsion(ff, mol, currentTorsionStamp);
120 mol->addTorsion(torsion);
121 }
122
123 //create inversions
124 Inversion* inversion;
125 InversionStamp* currentInversionStamp;
126 int nInversions = molStamp->getNInversions();
127 for (int i = 0; i < nInversions; ++i) {
128 currentInversionStamp = molStamp->getInversionStamp(i);
129 inversion = createInversion(ff, mol, currentInversionStamp);
130 if (inversion != NULL ) {
131 mol->addInversion(inversion);
132 }
133 }
134
135 //create cutoffGroups
136 CutoffGroup* cutoffGroup;
137 CutoffGroupStamp* currentCutoffGroupStamp;
138 int nCutoffGroups = molStamp->getNCutoffGroups();
139 for (int i = 0; i < nCutoffGroups; ++i) {
140 currentCutoffGroupStamp = molStamp->getCutoffGroupStamp(i);
141 cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp, localIndexMan);
142 mol->addCutoffGroup(cutoffGroup);
143 }
144
145 //every free atom is a cutoff group
146 std::vector<Atom*> freeAtoms;
147 std::vector<Atom*>::iterator ai;
148 std::vector<Atom*>::iterator fai;
149
150 //add all atoms into allAtoms set
151 for(atom = mol->beginAtom(fai); atom != NULL; atom = mol->nextAtom(fai)) {
152 freeAtoms.push_back(atom);
153 }
154
155 Molecule::CutoffGroupIterator ci;
156 CutoffGroup* cg;
157
158 for (cg = mol->beginCutoffGroup(ci); cg != NULL;
159 cg = mol->nextCutoffGroup(ci)) {
160
161 for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
162 //erase the atoms belong to cutoff groups from freeAtoms vector
163 freeAtoms.erase(std::remove(freeAtoms.begin(), freeAtoms.end(), atom),
164 freeAtoms.end());
165 }
166 }
167
168 // loop over the free atoms and then create one cutoff group for
169 // every single free atom
170
171 for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) {
172 cutoffGroup = createCutoffGroup(mol, *fai, localIndexMan);
173 mol->addCutoffGroup(cutoffGroup);
174 }
175 //create constraints
176 createConstraintPair(mol);
177 createConstraintElem(mol);
178
179 //the construction of this molecule is finished
180 mol->complete();
181
182 return mol;
183 }
184
185
186 Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol,
187 AtomStamp* stamp,
188 LocalIndexManager* localIndexMan) {
189 AtomType * atomType;
190 Atom* atom;
191
192 atomType = ff->getAtomType(stamp->getType());
193
194 if (atomType == NULL) {
195 sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
196 stamp->getType().c_str());
197
198 painCave.isFatal = 1;
199 simError();
200 }
201
202 //below code still have some kind of hard-coding smell
203 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");
209
210 painCave.isFatal = 1;
211 simError();
212 }
213
214 DirectionalAtom* dAtom;
215 dAtom = new DirectionalAtom(dAtomType);
216 atom = dAtom;
217 }
218 else{
219 atom = new Atom(atomType);
220 }
221
222 atom->setLocalIndex(localIndexMan->getNextAtomIndex());
223
224 return atom;
225 }
226
227 RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp,
228 Molecule* mol,
229 RigidBodyStamp* rbStamp,
230 LocalIndexManager* localIndexMan) {
231 Atom* atom;
232 int nAtoms;
233 Vector3d refCoor;
234 AtomStamp* atomStamp;
235
236 RigidBody* rb = new RigidBody();
237 nAtoms = rbStamp->getNMembers();
238 for (int i = 0; i < nAtoms; ++i) {
239 //rbStamp->getMember(i) return the local index of current atom
240 //inside the molecule. It is not the same as local index of
241 //atom which is the index of atom at DataStorage class
242 atom = mol->getAtomAt(rbStamp->getMemberAt(i));
243 atomStamp= molStamp->getAtomStamp(rbStamp->getMemberAt(i));
244 rb->addAtom(atom, atomStamp);
245 }
246
247 //after all of the atoms are added, we need to calculate the
248 //reference coordinates
249 rb->calcRefCoords();
250
251 //set the local index of this rigid body, global index will be set later
252 rb->setLocalIndex(localIndexMan->getNextRigidBodyIndex());
253
254 //the rule for naming rigidbody MoleculeName_RB_Integer
255 //The first part is the name of the molecule
256 //The second part is alway fixed as "RB"
257 //The third part is the index of the rigidbody defined in meta-data file
258 //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 = OpenMD_itoa(mol->getNRigidBodies(), 10);
261 rb->setType(mol->getType() + "_RB_" + s.c_str());
262
263 return rb;
264 }
265
266 Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol,
267 BondStamp* stamp) {
268 BondType* bondType;
269 Atom* atomA;
270 Atom* atomB;
271
272 atomA = mol->getAtomAt(stamp->getA());
273 atomB = mol->getAtomAt(stamp->getB());
274
275 assert( atomA && atomB);
276
277 bondType = ff->getBondType(atomA->getType(), atomB->getType());
278
279 if (bondType == NULL) {
280 sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
281 atomA->getType().c_str(),
282 atomB->getType().c_str());
283
284 painCave.isFatal = 1;
285 simError();
286 }
287 return new Bond(atomA, atomB, bondType);
288 }
289
290 Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol,
291 BendStamp* stamp) {
292 Bend* bend = NULL;
293 std::vector<int> bendAtoms = stamp->getMembers();
294 if (bendAtoms.size() == 3) {
295 Atom* atomA = mol->getAtomAt(bendAtoms[0]);
296 Atom* atomB = mol->getAtomAt(bendAtoms[1]);
297 Atom* atomC = mol->getAtomAt(bendAtoms[2]);
298
299 assert( atomA && atomB && atomC);
300
301 BendType* bendType = ff->getBendType(atomA->getType().c_str(),
302 atomB->getType().c_str(),
303 atomC->getType().c_str());
304
305 if (bendType == NULL) {
306 sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
307 atomA->getType().c_str(),
308 atomB->getType().c_str(),
309 atomC->getType().c_str());
310
311 painCave.isFatal = 1;
312 simError();
313 }
314
315 bend = new Bend(atomA, atomB, atomC, bendType);
316 } else if ( bendAtoms.size() == 2 && stamp->haveGhostVectorSource()) {
317 int ghostIndex = stamp->getGhostVectorSource();
318 int normalIndex = ghostIndex != bendAtoms[0] ? bendAtoms[0] : bendAtoms[1];
319 Atom* normalAtom = mol->getAtomAt(normalIndex) ;
320 DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
321 if (ghostAtom == NULL) {
322 sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
323 painCave.isFatal = 1;
324 simError();
325 }
326
327 BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
328
329 if (bendType == NULL) {
330 sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
331 normalAtom->getType().c_str(),
332 ghostAtom->getType().c_str(),
333 "GHOST");
334
335 painCave.isFatal = 1;
336 simError();
337 }
338
339 bend = new GhostBend(normalAtom, ghostAtom, bendType);
340
341 }
342
343 return bend;
344 }
345
346 Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol,
347 TorsionStamp* stamp) {
348
349 Torsion* torsion = NULL;
350 std::vector<int> torsionAtoms = stamp->getMembers();
351 if (torsionAtoms.size() < 3) {
352 return torsion;
353 }
354
355 Atom* atomA = mol->getAtomAt(torsionAtoms[0]);
356 Atom* atomB = mol->getAtomAt(torsionAtoms[1]);
357 Atom* atomC = mol->getAtomAt(torsionAtoms[2]);
358
359 if (torsionAtoms.size() == 4) {
360 Atom* atomD = mol->getAtomAt(torsionAtoms[3]);
361
362 assert(atomA && atomB && atomC && atomD);
363
364 TorsionType* torsionType = ff->getTorsionType(atomA->getType(),
365 atomB->getType(),
366 atomC->getType(),
367 atomD->getType());
368 if (torsionType == NULL) {
369 sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
370 atomA->getType().c_str(),
371 atomB->getType().c_str(),
372 atomC->getType().c_str(),
373 atomD->getType().c_str());
374
375 painCave.isFatal = 1;
376 simError();
377 }
378
379 torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);
380 }
381 else {
382
383 DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(stamp->getGhostVectorSource()));
384 if (dAtom == NULL) {
385 sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
386 painCave.isFatal = 1;
387 simError();
388 }
389
390 TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
391 atomC->getType(), "GHOST");
392
393 if (torsionType == NULL) {
394 sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
395 atomA->getType().c_str(),
396 atomB->getType().c_str(),
397 atomC->getType().c_str(),
398 "GHOST");
399
400 painCave.isFatal = 1;
401 simError();
402 }
403
404 torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);
405 }
406
407 return torsion;
408 }
409
410 Inversion* MoleculeCreator::createInversion(ForceField* ff, Molecule* mol,
411 InversionStamp* stamp) {
412
413 Inversion* inversion = NULL;
414 int center = stamp->getCenter();
415 std::vector<int> satellites = stamp->getSatellites();
416 if (satellites.size() != 3) {
417 return inversion;
418 }
419
420 Atom* atomA = mol->getAtomAt(center);
421 Atom* atomB = mol->getAtomAt(satellites[0]);
422 Atom* atomC = mol->getAtomAt(satellites[1]);
423 Atom* atomD = mol->getAtomAt(satellites[2]);
424
425 assert(atomA && atomB && atomC && atomD);
426
427 InversionType* inversionType = ff->getInversionType(atomA->getType(),
428 atomB->getType(),
429 atomC->getType(),
430 atomD->getType());
431
432 if (inversionType == NULL) {
433 sprintf(painCave.errMsg, "No Matching Inversion Type for[%s, %s, %s, %s]\n"
434 "\t(May not be a problem: not all inversions are parametrized)\n",
435 atomA->getType().c_str(),
436 atomB->getType().c_str(),
437 atomC->getType().c_str(),
438 atomD->getType().c_str());
439
440 painCave.isFatal = 0;
441 painCave.severity = OPENMD_INFO;
442 simError();
443 return NULL;
444 } else {
445
446 inversion = new Inversion(atomA, atomB, atomC, atomD, inversionType);
447 return inversion;
448 }
449 }
450
451
452 CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol,
453 CutoffGroupStamp* stamp,
454 LocalIndexManager* localIndexMan) {
455 int nAtoms;
456 CutoffGroup* cg;
457 Atom* atom;
458 cg = new CutoffGroup();
459
460 nAtoms = stamp->getNMembers();
461 for (int i =0; i < nAtoms; ++i) {
462 atom = mol->getAtomAt(stamp->getMemberAt(i));
463 assert(atom);
464 cg->addAtom(atom);
465 }
466
467 //set the local index of this cutoffGroup, global index will be set later
468 cg->setLocalIndex(localIndexMan->getNextCutoffGroupIndex());
469
470 return cg;
471 }
472
473 CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom,
474 LocalIndexManager* localIndexMan) {
475 CutoffGroup* cg;
476 cg = new CutoffGroup();
477 cg->addAtom(atom);
478
479 //set the local index of this cutoffGroup, global index will be set later
480 cg->setLocalIndex(localIndexMan->getNextCutoffGroupIndex());
481
482 return cg;
483 }
484
485 void MoleculeCreator::createConstraintPair(Molecule* mol) {
486
487 //add bond constraints
488 Molecule::BondIterator bi;
489 Bond* bond;
490 for (bond = mol->beginBond(bi); bond != NULL; bond = mol->nextBond(bi)) {
491
492 BondType* bt = bond->getBondType();
493
494 //class Parent1 {};
495 //class Child1 : public Parent {};
496 //class Child2 : public Parent {};
497 //Child1* ch1 = new Child1();
498 //Child2* ch2 = dynamic_cast<Child2*>(ch1);
499 //the dynamic_cast is succeed in above line. A compiler bug?
500
501 if (typeid(FixedBondType) == typeid(*bt)) {
502 FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt);
503
504 ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA());
505 ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB());
506 ConstraintPair* consPair = new ConstraintPair(consElemA, consElemB, fbt->getEquilibriumBondLength());
507 mol->addConstraintPair(consPair);
508 }
509 }
510
511 //rigidbody -- rigidbody constraint is not support yet
512 }
513
514 void MoleculeCreator::createConstraintElem(Molecule* mol) {
515
516 ConstraintPair* consPair;
517 Molecule::ConstraintPairIterator cpi;
518 std::set<StuntDouble*> sdSet;
519 for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
520
521 StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble();
522 if (sdSet.find(sdA) == sdSet.end()){
523 sdSet.insert(sdA);
524 mol->addConstraintElem(new ConstraintElem(sdA));
525 }
526
527 StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();
528 if (sdSet.find(sdB) == sdSet.end()){
529 sdSet.insert(sdB);
530 mol->addConstraintElem(new ConstraintElem(sdB));
531 }
532
533 }
534
535 }
536
537 }

Properties

Name Value
svn:executable *
svn:keywords Author Id Revision Date